options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

regression

normal linear model

ex5-1.stan

normal regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~normal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=normal_rng(m1[i],s);
  }
}
n=30
mdl=cmdstan_model('./ex5-1.stan')

single regression

tb=tibble(x=runif(n,0,9),
          y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)

qplot(data=tb,x,y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -8002.76 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -12.9393   0.000249625   0.000698878           1           1       52    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##    lp__     -12.94
##    b[1]       0.05
##    b[2]       0.99
##    s          0.93
##    y1[1]      2.68
##    y1[2]      3.09
##    y1[3]      3.94
##    y1[4]      5.56
##    y1[5]      5.71
##    y1[6]      5.54
##    y1[7]      4.13
##    y1[8]      3.11
##    y1[9]      0.37
##    y1[10]     8.16
##    y1[11]     6.22
##    y1[12]     4.65
##    y1[13]     0.93
##    y1[14]     2.82
##    y1[15]     2.32
##    y1[16]     3.06
##    y1[17]     6.35
##    y1[18]     9.67
##    y1[19]     9.02
##    y1[20]     3.89
##    y1[21]     6.47
##    y1[22]     5.50
##    y1[23]     0.79
##    y1[24]     2.50
##    y1[25]     5.07
##    y1[26]     0.58
##    y1[27]     9.13
##    y1[28]     0.87
##    y1[29]     2.85
##    y1[30]     2.76
##    m1[1]      0.68
##    m1[2]      2.84
##    m1[3]      4.15
##    m1[4]      5.05
##    m1[5]      6.02
##    m1[6]      5.07
##    m1[7]      3.79
##    m1[8]      0.89
##    m1[9]      1.51
##    m1[10]     8.86
##    m1[11]     5.30
##    m1[12]     4.44
##    m1[13]     0.95
##    m1[14]     3.59
##    m1[15]     2.09
##    m1[16]     2.28
##    m1[17]     7.40
##    m1[18]     8.43
##    m1[19]     8.64
##    m1[20]     2.88
##    m1[21]     6.00
##    m1[22]     5.47
##    m1[23]     0.96
##    m1[24]     2.81
##    m1[25]     5.82
##    m1[26]     1.49
##    m1[27]     8.31
##    m1[28]     1.27
##    m1[29]     1.41
##    m1[30]     4.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.62 -14.28 1.27 1.06 -17.12 -13.20 1.01      849     1150
##    b[1]     0.07   0.08 0.36 0.36  -0.53   0.65 1.01      396      545
##    b[2]     0.99   0.99 0.07 0.08   0.87   1.11 1.01      459      778
##    s        1.02   1.00 0.15 0.14   0.81   1.28 1.00     1039      948
##    y1[1]    0.68   0.65 1.04 1.03  -0.92   2.39 1.00     1794     1922
##    y1[2]    2.86   2.86 1.04 1.04   1.16   4.53 1.00     1981     1686
##    y1[3]    4.17   4.18 1.04 1.04   2.49   5.84 1.00     1965     1794
##    y1[4]    5.06   5.08 1.02 1.00   3.41   6.78 1.00     1947     1933
##    y1[5]    6.06   6.07 1.06 1.05   4.37   7.76 1.00     1780     1755
##    y1[6]    5.10   5.11 1.06 1.03   3.36   6.85 1.00     1979     1855
##    y1[7]    3.79   3.79 1.07 1.01   2.08   5.63 1.00     1831     1931
##    y1[8]    0.92   0.94 1.06 1.03  -0.88   2.69 1.00     1515     1848
##    y1[9]    1.51   1.49 1.03 1.01  -0.19   3.22 1.00     1686     1593
##    y1[10]   8.90   8.88 1.12 1.06   7.08  10.71 1.00     2112     1759
##    y1[11]   5.31   5.32 1.03 1.02   3.62   7.02 1.00     1789     1864
##    y1[12]   4.45   4.47 1.02 0.94   2.79   6.16 1.00     1832     1601
##    y1[13]   0.96   0.96 1.09 1.06  -0.82   2.73 1.00     1639     1775
##    y1[14]   3.62   3.62 1.05 1.05   1.90   5.35 1.00     1919     1846
##    y1[15]   2.06   2.05 1.05 0.97   0.28   3.80 1.00     1746     1879
##    y1[16]   2.32   2.33 1.05 1.01   0.60   4.05 1.00     2088     1981
##    y1[17]   7.41   7.38 1.08 1.06   5.64   9.18 1.00     1989     1971
##    y1[18]   8.43   8.40 1.09 1.07   6.72  10.18 1.00     2102     1917
##    y1[19]   8.65   8.64 1.11 1.09   6.91  10.51 1.00     1850     1863
##    y1[20]   2.87   2.90 1.04 0.99   1.15   4.55 1.00     1928     1890
##    y1[21]   5.97   5.96 1.06 1.01   4.21   7.72 1.00     1954     1859
##    y1[22]   5.54   5.55 1.05 0.98   3.81   7.23 1.00     2184     1957
##    y1[23]   0.96   0.99 1.08 1.05  -0.83   2.72 1.00     1808     1717
##    y1[24]   2.81   2.83 1.07 0.99   1.07   4.61 1.00     1748     1984
##    y1[25]   5.82   5.82 1.06 1.02   4.13   7.54 1.00     1925     1742
##    y1[26]   1.50   1.49 1.05 1.04  -0.25   3.21 1.00     1397     1597
##    y1[27]   8.31   8.29 1.09 1.07   6.49  10.08 1.00     1848     2104
##    y1[28]   1.25   1.25 1.08 1.05  -0.58   3.02 1.00     1815     1931
##    y1[29]   1.42   1.44 1.07 1.04  -0.39   3.18 1.00     1660     1913
##    y1[30]   4.00   4.01 1.08 1.08   2.19   5.74 1.00     1836     1903
##    m1[1]    0.70   0.71 0.32 0.32   0.16   1.21 1.01      410      474
##    m1[2]    2.85   2.85 0.21 0.21   2.50   3.20 1.01      632      983
##    m1[3]    4.16   4.15 0.19 0.19   3.85   4.47 1.00     1139     1265
##    m1[4]    5.05   5.05 0.20 0.20   4.73   5.39 1.00     1730     1376
##    m1[5]    6.02   6.01 0.24 0.23   5.65   6.40 1.00     1711     1540
##    m1[6]    5.08   5.07 0.20 0.20   4.75   5.41 1.00     1738     1376
##    m1[7]    3.80   3.80 0.19 0.19   3.49   4.12 1.00      940     1286
##    m1[8]    0.91   0.91 0.31 0.31   0.39   1.39 1.01      417      477
##    m1[9]    1.53   1.54 0.27 0.27   1.08   1.97 1.01      449      540
##    m1[10]   8.86   8.86 0.40 0.40   8.20   9.51 1.01      798     1184
##    m1[11]   5.30   5.29 0.21 0.20   4.96   5.65 1.00     1794     1400
##    m1[12]   4.44   4.44 0.19 0.19   4.13   4.76 1.00     1335     1302
##    m1[13]   0.96   0.97 0.30 0.30   0.45   1.45 1.01      420      477
##    m1[14]   3.60   3.59 0.19 0.19   3.28   3.92 1.00      855     1249
##    m1[15]   2.10   2.11 0.24 0.24   1.70   2.50 1.01      498      639
##    m1[16]   2.29   2.29 0.23 0.24   1.90   2.67 1.01      522      726
##    m1[17]   7.40   7.40 0.31 0.31   6.90   7.90 1.00     1032     1393
##    m1[18]   8.43   8.43 0.37 0.37   7.82   9.03 1.01      847     1180
##    m1[19]   8.64   8.64 0.39 0.39   8.00   9.26 1.01      822     1166
##    m1[20]   2.89   2.89 0.21 0.21   2.54   3.24 1.01      640      982
##    m1[21]   6.00   6.00 0.24 0.23   5.63   6.38 1.00     1717     1540
##    m1[22]   5.48   5.47 0.21 0.21   5.14   5.83 1.00     1823     1421
##    m1[23]   0.97   0.98 0.30 0.30   0.47   1.45 1.01      421      455
##    m1[24]   2.82   2.81 0.21 0.21   2.46   3.17 1.01      622      907
##    m1[25]   5.82   5.81 0.23 0.22   5.46   6.20 1.00     1778     1502
##    m1[26]   1.50   1.51 0.27 0.27   1.05   1.94 1.01      448      540
##    m1[27]   8.31   8.30 0.36 0.36   7.71   8.90 1.01      865     1154
##    m1[28]   1.29   1.30 0.28 0.28   0.82   1.75 1.01      434      539
##    m1[29]   1.42   1.43 0.28 0.27   0.96   1.87 1.01      441      507
##    m1[30]   4.01   4.00 0.19 0.19   3.70   4.32 1.00     1050     1266

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)

par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

multiple regression

tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -77.1658 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26       -13.335   0.000299111    0.00224091           1           1       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.34
##    b[1]       0.27
##    b[2]       1.04
##    b[3]      -1.06
##    s          0.95
##    y1[1]      4.97
##    y1[2]     -2.17
##    y1[3]     -0.15
##    y1[4]      0.72
##    y1[5]     -6.24
##    y1[6]      0.29
##    y1[7]      2.00
##    y1[8]     -2.90
##    y1[9]     -6.79
##    y1[10]    -0.56
##    y1[11]    -1.64
##    y1[12]     1.28
##    y1[13]     5.10
##    y1[14]     2.53
##    y1[15]     0.93
##    y1[16]     0.31
##    y1[17]     1.55
##    y1[18]     4.51
##    y1[19]    -3.39
##    y1[20]     1.61
##    y1[21]     3.69
##    y1[22]     5.02
##    y1[23]    -1.66
##    y1[24]     1.26
##    y1[25]    -1.39
##    y1[26]    -5.03
##    y1[27]    -1.85
##    y1[28]     5.50
##    y1[29]     7.91
##    y1[30]     2.80
##    m1[1]      5.18
##    m1[2]     -2.42
##    m1[3]     -0.93
##    m1[4]     -0.14
##    m1[5]     -6.83
##    m1[6]      1.33
##    m1[7]      2.48
##    m1[8]     -2.71
##    m1[9]     -6.60
##    m1[10]    -1.29
##    m1[11]    -0.59
##    m1[12]     1.12
##    m1[13]     3.26
##    m1[14]     1.60
##    m1[15]     0.39
##    m1[16]     1.86
##    m1[17]     0.32
##    m1[18]     4.74
##    m1[19]    -3.54
##    m1[20]    -1.26
##    m1[21]     3.51
##    m1[22]     4.28
##    m1[23]    -2.33
##    m1[24]     1.94
##    m1[25]    -1.03
##    m1[26]    -4.36
##    m1[27]    -2.96
##    m1[28]     6.75
##    m1[29]     7.23
##    m1[30]     3.18
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.58 -15.19 1.58 1.31 -18.60 -13.83 1.01      635      906
##    b[1]     0.28   0.29 0.41 0.42  -0.40   0.95 1.01      403      426
##    b[2]     1.04   1.04 0.07 0.07   0.92   1.16 1.00     1226     1018
##    b[3]    -1.06  -1.06 0.07 0.07  -1.18  -0.93 1.00      740     1105
##    s        1.05   1.03 0.16 0.15   0.82   1.34 1.00      941     1167
##    y1[1]    5.16   5.14 1.09 1.03   3.41   6.90 1.00     2029     1710
##    y1[2]   -2.40  -2.43 1.14 1.07  -4.24  -0.49 1.00     1895     1802
##    y1[3]   -0.94  -0.92 1.12 1.06  -2.75   0.92 1.00     1413     1880
##    y1[4]   -0.12  -0.13 1.09 1.09  -1.85   1.66 1.00     1711     1951
##    y1[5]   -6.81  -6.83 1.15 1.13  -8.59  -4.86 1.00     2099     1540
##    y1[6]    1.33   1.33 1.11 1.07  -0.53   3.13 1.00     1954     1548
##    y1[7]    2.47   2.47 1.13 1.09   0.63   4.32 1.00     1628     1900
##    y1[8]   -2.72  -2.73 1.10 1.08  -4.43  -0.95 1.00     2006     1785
##    y1[9]   -6.56  -6.55 1.12 1.09  -8.45  -4.72 1.00     2045     1934
##    y1[10]  -1.26  -1.23 1.08 1.06  -3.01   0.48 1.00     1782     1838
##    y1[11]  -0.57  -0.57 1.08 1.09  -2.38   1.18 1.00     1926     1737
##    y1[12]   1.07   1.06 1.13 1.11  -0.72   2.94 1.00     1541     1843
##    y1[13]   3.24   3.25 1.09 1.07   1.44   5.03 1.00     1910     1805
##    y1[14]   1.55   1.55 1.08 1.07  -0.19   3.36 1.00     1950     1968
##    y1[15]   0.38   0.37 1.08 1.05  -1.41   2.17 1.00     2062     1961
##    y1[16]   1.89   1.91 1.11 1.05   0.07   3.68 1.00     1537     1767
##    y1[17]   0.35   0.35 1.16 1.13  -1.52   2.24 1.00     1697     1786
##    y1[18]   4.71   4.74 1.10 1.06   2.93   6.46 1.00     1778     1778
##    y1[19]  -3.50  -3.52 1.10 1.12  -5.29  -1.70 1.00     1873     1840
##    y1[20]  -1.27  -1.30 1.12 1.10  -3.10   0.61 1.00     2058     2020
##    y1[21]   3.55   3.56 1.12 1.11   1.69   5.36 1.00     1792     1895
##    y1[22]   4.25   4.28 1.12 1.05   2.34   6.09 1.00     1865     1822
##    y1[23]  -2.33  -2.34 1.10 1.08  -4.12  -0.53 1.00     2027     1930
##    y1[24]   1.92   1.95 1.11 1.07   0.15   3.74 1.00     1736     1774
##    y1[25]  -1.08  -1.09 1.10 1.08  -2.87   0.73 1.00     1790     1875
##    y1[26]  -4.35  -4.35 1.13 1.09  -6.22  -2.54 1.00     1797     1857
##    y1[27]  -2.98  -2.97 1.09 1.07  -4.79  -1.24 1.00     2116     1916
##    y1[28]   6.72   6.70 1.12 1.08   4.97   8.59 1.00     2012     1881
##    y1[29]   7.23   7.19 1.15 1.12   5.42   9.17 1.00     1761     1684
##    y1[30]   3.17   3.18 1.10 1.05   1.41   5.04 1.00     1867     1871
##    m1[1]    5.17   5.17 0.34 0.34   4.62   5.72 1.00     1674     1279
##    m1[2]   -2.42  -2.42 0.37 0.35  -3.00  -1.81 1.01     1079     1195
##    m1[3]   -0.93  -0.91 0.35 0.35  -1.50  -0.37 1.01      441      451
##    m1[4]   -0.14  -0.13 0.29 0.29  -0.62   0.32 1.01      437      467
##    m1[5]   -6.83  -6.84 0.46 0.45  -7.56  -6.07 1.01     1993     1324
##    m1[6]    1.33   1.33 0.33 0.33   0.79   1.84 1.01     1256     1088
##    m1[7]    2.48   2.49 0.34 0.35   1.92   3.03 1.01      428      554
##    m1[8]   -2.70  -2.71 0.26 0.25  -3.12  -2.27 1.00     1478     1354
##    m1[9]   -6.59  -6.60 0.45 0.43  -7.31  -5.87 1.01     1997     1334
##    m1[10]  -1.28  -1.27 0.24 0.24  -1.70  -0.90 1.01      614      721
##    m1[11]  -0.59  -0.58 0.22 0.21  -0.95  -0.23 1.01     1970     1414
##    m1[12]   1.11   1.11 0.39 0.38   0.46   1.71 1.01     1007     1045
##    m1[13]   3.25   3.27 0.32 0.31   2.71   3.76 1.01     1735     1326
##    m1[14]   1.60   1.60 0.22 0.23   1.21   1.94 1.01     1723     1468
##    m1[15]   0.38   0.39 0.26 0.26  -0.04   0.78 1.01     1680     1298
##    m1[16]   1.86   1.87 0.36 0.37   1.26   2.45 1.01      406      459
##    m1[17]   0.32   0.32 0.43 0.41  -0.40   0.99 1.01      871      937
##    m1[18]   4.73   4.73 0.34 0.33   4.19   5.27 1.00     1924     1518
##    m1[19]  -3.53  -3.54 0.31 0.30  -4.05  -3.02 1.00      995     1248
##    m1[20]  -1.26  -1.26 0.22 0.21  -1.62  -0.90 1.01     1676     1222
##    m1[21]   3.50   3.51 0.36 0.35   2.90   4.06 1.01     1511     1267
##    m1[22]   4.28   4.28 0.32 0.31   3.74   4.81 1.00      652      967
##    m1[23]  -2.32  -2.32 0.26 0.25  -2.76  -1.89 1.01     2143     1391
##    m1[24]   1.94   1.94 0.35 0.35   1.36   2.51 1.01      410      490
##    m1[25]  -1.03  -1.01 0.36 0.36  -1.61  -0.45 1.01      442      425
##    m1[26]  -4.35  -4.35 0.36 0.36  -4.94  -3.76 1.00      914     1268
##    m1[27]  -2.96  -2.96 0.27 0.26  -3.40  -2.50 1.00     1998     1402
##    m1[28]   6.74   6.74 0.41 0.40   6.09   7.41 1.00     1605     1304
##    m1[29]   7.22   7.22 0.44 0.43   6.53   7.94 1.00     1098     1248
##    m1[30]   3.18   3.18 0.27 0.27   2.72   3.62 1.00      640      949

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANOVA

tb=tibble(c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)

qplot(data=tb,c,y,geom='boxplot')

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -41.02 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -15.6149   0.000624738   0.000527558           1           1       20    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -15.61
##    b[1]      -0.23
##    b[2]       2.36
##    b[3]      -1.82
##    s          1.02
##    y1[1]      2.84
##    y1[2]     -1.61
##    y1[3]     -1.08
##    y1[4]      1.85
##    y1[5]     -0.88
##    y1[6]      3.70
##    y1[7]     -0.40
##    y1[8]     -1.28
##    y1[9]      0.40
##    y1[10]     3.53
##    y1[11]    -0.22
##    y1[12]    -2.51
##    y1[13]    -0.41
##    y1[14]     1.97
##    y1[15]     2.60
##    y1[16]     2.48
##    y1[17]    -0.94
##    y1[18]     3.08
##    y1[19]    -2.31
##    y1[20]     0.47
##    y1[21]    -2.97
##    y1[22]    -2.25
##    y1[23]    -2.68
##    y1[24]    -1.02
##    y1[25]     1.34
##    y1[26]    -0.44
##    y1[27]     0.33
##    y1[28]     0.18
##    y1[29]    -1.39
##    y1[30]     2.63
##    m1[1]      2.13
##    m1[2]     -2.05
##    m1[3]     -0.23
##    m1[4]      2.13
##    m1[5]     -0.23
##    m1[6]      2.13
##    m1[7]     -0.23
##    m1[8]     -2.05
##    m1[9]     -0.23
##    m1[10]     2.13
##    m1[11]    -0.23
##    m1[12]    -2.05
##    m1[13]    -0.23
##    m1[14]     2.13
##    m1[15]     2.13
##    m1[16]     2.13
##    m1[17]    -0.23
##    m1[18]     2.13
##    m1[19]    -2.05
##    m1[20]    -0.23
##    m1[21]    -2.05
##    m1[22]    -2.05
##    m1[23]    -2.05
##    m1[24]    -2.05
##    m1[25]     2.13
##    m1[26]    -0.23
##    m1[27]    -0.23
##    m1[28]    -0.23
##    m1[29]    -2.05
##    m1[30]     2.13
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.78 -17.43 1.55 1.31 -20.81 -15.97 1.00      742     1057
##    b[1]    -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    b[2]     2.35   2.34 0.49 0.48   1.57   3.15 1.00      718      781
##    b[3]    -1.83  -1.83 0.55 0.52  -2.75  -0.92 1.00      635      740
##    s        1.13   1.11 0.17 0.16   0.90   1.42 1.00     1894     1741
##    y1[1]    2.11   2.14 1.21 1.17   0.15   4.07 1.00     1957     1893
##    y1[2]   -2.09  -2.07 1.20 1.16  -4.05  -0.17 1.00     1709     1833
##    y1[3]   -0.20  -0.18 1.21 1.18  -2.25   1.77 1.00     1740     1921
##    y1[4]    2.15   2.17 1.22 1.20   0.14   4.15 1.00     1675     1536
##    y1[5]   -0.19  -0.20 1.20 1.19  -2.11   1.81 1.00     1596     1697
##    y1[6]    2.12   2.11 1.21 1.22   0.14   4.17 1.00     1838     1999
##    y1[7]   -0.22  -0.18 1.17 1.16  -2.19   1.67 1.00     1862     1735
##    y1[8]   -2.07  -2.11 1.22 1.17  -4.07  -0.05 1.00     1883     1860
##    y1[9]   -0.21  -0.20 1.19 1.21  -2.19   1.76 1.00     1426     1901
##    y1[10]   2.16   2.16 1.23 1.20   0.16   4.11 1.00     1743     1811
##    y1[11]  -0.24  -0.24 1.20 1.12  -2.15   1.76 1.00     1508     1956
##    y1[12]  -2.06  -2.04 1.21 1.19  -4.08  -0.10 1.00     1892     1700
##    y1[13]  -0.25  -0.25 1.18 1.12  -2.17   1.75 1.00     1747     1973
##    y1[14]   2.11   2.13 1.18 1.17   0.15   4.01 1.00     2076     1973
##    y1[15]   2.12   2.10 1.17 1.13   0.20   4.02 1.00     2012     1994
##    y1[16]   2.15   2.15 1.18 1.14   0.22   4.13 1.00     1913     1643
##    y1[17]  -0.23  -0.22 1.19 1.18  -2.22   1.71 1.00     1751     1915
##    y1[18]   2.14   2.11 1.24 1.26   0.12   4.19 1.00     1924     1855
##    y1[19]  -2.05  -2.04 1.19 1.16  -4.01  -0.11 1.00     1835     1806
##    y1[20]  -0.21  -0.21 1.19 1.17  -2.12   1.79 1.00     1812     1774
##    y1[21]  -2.06  -2.05 1.24 1.18  -4.07   0.00 1.00     1849     1744
##    y1[22]  -2.08  -2.10 1.21 1.13  -4.07  -0.09 1.00     1695     1938
##    y1[23]  -2.08  -2.06 1.22 1.18  -4.08  -0.07 1.00     1676     1680
##    y1[24]  -2.05  -2.03 1.20 1.19  -4.05  -0.05 1.00     1855     2013
##    y1[25]   2.12   2.11 1.21 1.23   0.15   4.10 1.00     1975     1931
##    y1[26]  -0.22  -0.18 1.19 1.11  -2.26   1.70 1.00     1784     1899
##    y1[27]  -0.24  -0.21 1.19 1.18  -2.23   1.68 1.00     1864     1846
##    y1[28]  -0.25  -0.24 1.21 1.18  -2.22   1.77 1.00     1737     1910
##    y1[29]  -2.02  -2.02 1.20 1.19  -3.96  -0.03 1.00     1832     1971
##    y1[30]   2.12   2.08 1.21 1.15   0.15   4.22 1.00     1867     1963
##    m1[1]    2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[2]   -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[3]   -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[4]    2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[5]   -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[6]    2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[7]   -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[8]   -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[9]   -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[10]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[11]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[12]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[13]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[14]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[15]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[16]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[17]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[18]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[19]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[20]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[21]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[22]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[23]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[24]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[25]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387
##    m1[26]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[27]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[28]  -0.22  -0.21 0.35 0.33  -0.79   0.36 1.00      709      916
##    m1[29]  -2.05  -2.05 0.40 0.38  -2.71  -1.41 1.00     1205     1120
##    m1[30]   2.13   2.13 0.36 0.35   1.53   2.70 1.00     1404     1387

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANCOVA

tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))

f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)

lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])

qplot(data=tb,x,y,shape=c,size=I(2))

plot(tb$x,tb$y,col=1+lv[factor(tb$c)])

qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -5073.62 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       38      -14.9114   4.98066e-05    0.00237982      0.2792           1       69    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -14.91
##    b[1]      -0.25
##    b[2]       1.03
##    b[3]       2.00
##    b[4]      -2.81
##    s          1.00
##    y1[1]      8.11
##    y1[2]     -0.56
##    y1[3]      1.60
##    y1[4]      4.14
##    y1[5]      8.45
##    y1[6]      6.87
##    y1[7]      4.97
##    y1[8]      2.45
##    y1[9]      1.75
##    y1[10]     8.35
##    y1[11]     4.08
##    y1[12]     4.12
##    y1[13]     5.94
##    y1[14]    -2.97
##    y1[15]     6.52
##    y1[16]     1.97
##    y1[17]    10.81
##    y1[18]     6.85
##    y1[19]     1.38
##    y1[20]     2.96
##    y1[21]     4.67
##    y1[22]     5.19
##    y1[23]    12.04
##    y1[24]     4.15
##    y1[25]     1.46
##    y1[26]     5.92
##    y1[27]     2.57
##    y1[28]     3.58
##    y1[29]     1.33
##    y1[30]     3.03
##    m1[1]      9.23
##    m1[2]     -2.68
##    m1[3]      1.67
##    m1[4]      4.12
##    m1[5]      8.70
##    m1[6]      8.27
##    m1[7]      5.42
##    m1[8]      1.85
##    m1[9]      1.85
##    m1[10]     9.26
##    m1[11]     2.42
##    m1[12]     4.36
##    m1[13]     7.03
##    m1[14]    -2.14
##    m1[15]     6.39
##    m1[16]     1.65
##    m1[17]     8.92
##    m1[18]     7.32
##    m1[19]     1.90
##    m1[20]     1.51
##    m1[21]     6.72
##    m1[22]     5.22
##    m1[23]    10.57
##    m1[24]     3.17
##    m1[25]     0.32
##    m1[26]     5.47
##    m1[27]     2.65
##    m1[28]     5.07
##    m1[29]     1.26
##    m1[30]     2.92
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.69 -17.31 1.76 1.55 -20.89 -15.55 1.00      609     1054
##    b[1]    -0.25  -0.27 0.46 0.43  -0.98   0.53 1.00      726      988
##    b[2]     1.03   1.03 0.08 0.08   0.90   1.16 1.00     1598     1659
##    b[3]     2.00   2.02 0.49 0.47   1.18   2.82 1.00      872      906
##    b[4]    -2.83  -2.82 0.60 0.57  -3.82  -1.81 1.01      471      821
##    s        1.12   1.11 0.17 0.16   0.89   1.42 1.00     1720     1177
##    y1[1]    9.25   9.25 1.18 1.19   7.31  11.18 1.00     2239     1904
##    y1[2]   -2.64  -2.63 1.27 1.23  -4.69  -0.57 1.00     1880     1930
##    y1[3]    1.67   1.66 1.19 1.19  -0.25   3.61 1.00     1726     2013
##    y1[4]    4.12   4.10 1.20 1.19   2.17   6.07 1.00     2015     1931
##    y1[5]    8.69   8.71 1.20 1.12   6.65  10.62 1.00     1946     1918
##    y1[6]    8.25   8.26 1.17 1.09   6.32  10.13 1.00     1852     2014
##    y1[7]    5.43   5.42 1.26 1.19   3.38   7.46 1.00     1598     1878
##    y1[8]    1.86   1.84 1.26 1.24  -0.19   3.97 1.00     1918     1979
##    y1[9]    1.81   1.82 1.19 1.10  -0.12   3.75 1.00     1832     1631
##    y1[10]   9.24   9.24 1.19 1.21   7.31  11.15 1.00     1994     1881
##    y1[11]   2.45   2.46 1.25 1.20   0.37   4.53 1.00     1876     1854
##    y1[12]   4.30   4.30 1.17 1.13   2.37   6.24 1.00     1795     2015
##    y1[13]   7.03   7.02 1.16 1.13   5.16   8.97 1.00     1875     1974
##    y1[14]  -2.18  -2.19 1.22 1.22  -4.18  -0.16 1.00     1796     2031
##    y1[15]   6.38   6.39 1.22 1.16   4.31   8.39 1.00     1818     1685
##    y1[16]   1.68   1.70 1.22 1.18  -0.35   3.65 1.00     1813     1937
##    y1[17]   8.93   8.94 1.19 1.13   7.03  10.90 1.00     1608     1743
##    y1[18]   7.30   7.30 1.14 1.08   5.42   9.24 1.00     1860     1927
##    y1[19]   1.90   1.87 1.21 1.19   0.01   3.87 1.00     1957     1684
##    y1[20]   1.49   1.51 1.22 1.21  -0.52   3.43 1.00     1685     1932
##    y1[21]   6.74   6.73 1.24 1.17   4.72   8.85 1.00     1863     1933
##    y1[22]   5.21   5.22 1.22 1.17   3.14   7.21 1.00     1712     1540
##    y1[23]  10.59  10.61 1.23 1.15   8.52  12.57 1.00     1669     1582
##    y1[24]   3.18   3.19 1.23 1.18   1.06   5.15 1.00     1960     1932
##    y1[25]   0.31   0.32 1.25 1.20  -1.72   2.33 1.00     1531     1836
##    y1[26]   5.49   5.49 1.16 1.17   3.63   7.43 1.00     2064     2000
##    y1[27]   2.65   2.64 1.16 1.11   0.76   4.58 1.00     1795     1951
##    y1[28]   5.03   5.04 1.18 1.14   3.07   6.98 1.00     2029     1896
##    y1[29]   1.24   1.27 1.18 1.18  -0.72   3.22 1.00     1753     1912
##    y1[30]   2.93   2.94 1.21 1.17   0.94   4.87 1.00     2043     1835
##    m1[1]    9.24   9.24 0.38 0.37   8.62   9.85 1.00     1620     1447
##    m1[2]   -2.70  -2.70 0.53 0.50  -3.54  -1.82 1.00      928     1062
##    m1[3]    1.67   1.64 0.39 0.37   1.04   2.30 1.00      643      877
##    m1[4]    4.12   4.12 0.38 0.38   3.49   4.74 1.00     1435     1079
##    m1[5]    8.70   8.70 0.36 0.35   8.11   9.27 1.00     1597     1470
##    m1[6]    8.27   8.27 0.34 0.33   7.71   8.83 1.00     1578     1469
##    m1[7]    5.41   5.42 0.61 0.58   4.42   6.39 1.00     1051     1465
##    m1[8]    1.84   1.85 0.50 0.50   1.04   2.65 1.00     1451      881
##    m1[9]    1.85   1.83 0.38 0.36   1.23   2.48 1.00      632      866
##    m1[10]   9.26   9.27 0.38 0.37   8.64   9.88 1.00     1622     1447
##    m1[11]   2.40   2.41 0.49 0.49   1.61   3.18 1.01      883     1215
##    m1[12]   4.36   4.35 0.38 0.37   3.74   4.99 1.00      705      851
##    m1[13]   7.03   7.03 0.32 0.31   6.52   7.56 1.00     1511     1567
##    m1[14]  -2.15  -2.15 0.51 0.50  -2.97  -1.31 1.00      897      988
##    m1[15]   6.40   6.40 0.44 0.42   5.65   7.11 1.00     1050     1270
##    m1[16]   1.64   1.65 0.48 0.47   0.86   2.41 1.01      858     1251
##    m1[17]   8.92   8.93 0.36 0.36   8.32   9.51 1.00     1607     1404
##    m1[18]   7.32   7.32 0.32 0.31   6.79   7.84 1.00     1527     1541
##    m1[19]   1.90   1.88 0.38 0.36   1.28   2.52 1.00      631      806
##    m1[20]   1.50   1.51 0.48 0.47   0.72   2.26 1.01      854     1348
##    m1[21]   6.72   6.73 0.46 0.44   5.95   7.45 1.00     1090     1263
##    m1[22]   5.22   5.22 0.40 0.38   4.57   5.89 1.00      786     1087
##    m1[23]  10.58  10.58 0.44 0.44   9.84  11.29 1.00     1655     1516
##    m1[24]   3.17   3.17 0.43 0.43   2.47   3.86 1.00     1434      961
##    m1[25]   0.31   0.29 0.43 0.41  -0.38   1.05 1.00      705      958
##    m1[26]   5.47   5.47 0.33 0.33   4.90   6.01 1.00     1455     1212
##    m1[27]   2.65   2.63 0.37 0.36   2.05   3.27 1.00      626      882
##    m1[28]   5.07   5.07 0.35 0.34   4.49   5.63 1.00     1445     1195
##    m1[29]   1.26   1.24 0.40 0.38   0.62   1.93 1.00      666      889
##    m1[30]   2.92   2.93 0.44 0.44   2.20   3.63 1.00     1436      966

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
  geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

interaction of variable

n=50
mdl=cmdstan_model('./ex5-1.stan')

tb=tibble(x=runif(n,-3,3),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,x,y,col=ca),
             qplot(data=tb,x,y,col=cb),ncol=2)

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)

fn(f0)
## Initial log joint probability = -2013.08 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13      -20.4747   7.45039e-05   0.000373311      0.8734      0.8734       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.59 -22.30 1.43 1.26 -25.35 -20.94 1.00      811     1095
##    b[1]     0.33   0.33 0.17 0.16   0.04   0.61 1.00      765      972
##    b[2]     1.22   1.22 0.08 0.08   1.09   1.36 1.00     2614     1715
##    b[3]     0.44   0.44 0.29 0.27  -0.03   0.93 1.00      630      843
##    s        0.96   0.96 0.10 0.10   0.81   1.15 1.00     1850     1570
##    y1[1]    1.31   1.27 0.97 0.97  -0.29   2.90 1.00     2035     2057
##    y1[2]    1.07   1.04 0.97 0.94  -0.49   2.68 1.00     2021     1967
##    y1[3]    3.16   3.16 0.99 0.95   1.54   4.75 1.00     1971     1854
##    y1[4]    2.68   2.72 1.01 0.99   0.99   4.38 1.00     2033     1988
##    y1[5]   -2.98  -2.98 1.02 1.01  -4.60  -1.34 1.00     1845     1930
##    y1[6]   -1.62  -1.59 1.00 0.96  -3.32   0.05 1.00     1737     1912
##    y1[7]    0.58   0.59 0.96 0.93  -0.94   2.18 1.00     1863     1893
##    y1[8]   -1.34  -1.33 1.00 0.98  -3.05   0.30 1.00     1954     1848
##    y1[9]   -0.85  -0.86 0.99 0.98  -2.45   0.76 1.00     1889     1817
##    y1[10]   2.35   2.37 0.99 0.97   0.72   4.00 1.00     2050     1973
##    y1[11]   2.37   2.39 1.01 1.00   0.72   4.04 1.00     2154     1907
##    y1[12]   0.02   0.01 0.99 0.99  -1.60   1.65 1.00     1931     1971
##    y1[13]  -1.85  -1.87 1.02 1.04  -3.55  -0.14 1.00     2077     2055
##    y1[14]   3.21   3.19 1.01 0.97   1.61   4.87 1.00     2107     2048
##    y1[15]  -3.27  -3.26 1.01 0.99  -5.01  -1.57 1.00     1923     1660
##    y1[16]   0.44   0.45 0.97 0.94  -1.12   2.01 1.00     1888     1959
##    y1[17]   2.79   2.78 1.01 1.01   1.16   4.44 1.00     1876     1883
##    y1[18]  -1.56  -1.58 0.99 0.97  -3.15   0.07 1.00     1901     1880
##    y1[19]   3.33   3.35 1.01 1.02   1.67   5.00 1.00     1876     1822
##    y1[20]   3.65   3.65 1.00 0.99   2.00   5.29 1.00     1931     1844
##    y1[21]   2.17   2.17 0.98 0.99   0.57   3.77 1.00     1865     1959
##    y1[22]  -1.83  -1.85 1.03 0.99  -3.55  -0.13 1.00     1811     1933
##    y1[23]   1.58   1.60 0.98 0.95  -0.03   3.23 1.00     1898     2103
##    y1[24]   0.41   0.42 1.00 0.95  -1.24   2.02 1.00     1883     1904
##    y1[25]   1.38   1.39 0.95 0.97  -0.23   2.91 1.00     1757     1772
##    y1[26]  -0.20  -0.19 1.00 0.99  -1.88   1.47 1.00     1829     1977
##    y1[27]   2.10   2.14 0.98 0.97   0.43   3.68 1.00     1807     1970
##    y1[28]   0.96   0.93 0.97 0.99  -0.65   2.55 1.00     2021     2000
##    y1[29]  -2.17  -2.19 1.02 1.02  -3.82  -0.43 1.00     1709     1790
##    y1[30]   2.56   2.58 0.99 0.94   0.89   4.18 1.00     1719     1839
##    y1[31]   2.06   2.04 1.00 0.95   0.41   3.71 1.00     1821     1931
##    y1[32]   4.29   4.30 1.01 1.00   2.58   5.93 1.00     1864     1974
##    y1[33]   3.94   3.96 1.03 1.02   2.20   5.58 1.00     1816     1700
##    y1[34]   2.58   2.58 1.01 0.98   0.89   4.23 1.00     1460     2001
##    y1[35]   1.14   1.16 1.00 0.94  -0.52   2.75 1.00     1970     1558
##    y1[36]  -0.20  -0.18 1.00 0.99  -1.88   1.45 1.00     1934     2010
##    y1[37]  -0.61  -0.63 0.96 0.92  -2.16   1.00 1.00     1975     1968
##    y1[38]  -1.51  -1.51 0.99 0.95  -3.14   0.11 1.00     1928     2054
##    y1[39]   3.14   3.13 1.00 0.98   1.50   4.78 1.00     1807     1955
##    y1[40]  -2.97  -2.96 1.00 1.00  -4.56  -1.38 1.00     1854     1917
##    y1[41]   2.88   2.93 1.01 1.00   1.20   4.54 1.00     1864     1819
##    y1[42]   1.99   1.98 0.98 0.97   0.41   3.61 1.00     2006     1933
##    y1[43]   0.37   0.39 1.00 0.97  -1.26   1.95 1.00     1894     1968
##    y1[44]   3.92   3.93 1.02 1.02   2.17   5.54 1.00     2036     2011
##    y1[45]  -1.41  -1.42 1.01 0.99  -3.08   0.29 1.00     1837     1807
##    y1[46]   2.57   2.58 1.00 0.99   0.99   4.24 1.00     1807     1852
##    y1[47]  -0.48  -0.51 1.00 0.98  -2.14   1.13 1.00     1797     2057
##    y1[48]   1.64   1.64 0.98 0.92  -0.05   3.28 1.00     1663     1729
##    y1[49]   2.88   2.88 1.01 0.98   1.23   4.51 1.00     1860     1783
##    y1[50]  -0.24  -0.25 0.99 0.94  -1.82   1.36 1.00     1818     1777
##    m1[1]    1.27   1.28 0.17 0.17   0.99   1.57 1.00      817     1207
##    m1[2]    1.07   1.08 0.17 0.16   0.79   1.36 1.00      793     1068
##    m1[3]    3.17   3.17 0.26 0.25   2.73   3.59 1.00     1189     1457
##    m1[4]    2.69   2.69 0.25 0.23   2.29   3.09 1.00     1109     1431
##    m1[5]   -2.94  -2.94 0.30 0.29  -3.43  -2.45 1.00     1374     1492
##    m1[6]   -1.63  -1.62 0.29 0.29  -2.10  -1.16 1.00     1301     1471
##    m1[7]    0.62   0.62 0.17 0.16   0.34   0.90 1.00      762     1011
##    m1[8]   -1.36  -1.37 0.22 0.21  -1.74  -1.00 1.00     1039      989
##    m1[9]   -0.83  -0.83 0.20 0.19  -1.16  -0.50 1.00      930     1070
##    m1[10]   2.35   2.35 0.20 0.20   2.02   2.70 1.00     1061     1330
##    m1[11]   2.38   2.37 0.24 0.23   1.99   2.77 1.00     1062     1390
##    m1[12]   0.02   0.02 0.18 0.17  -0.28   0.31 1.00      786      986
##    m1[13]  -1.84  -1.83 0.30 0.30  -2.33  -1.36 1.00     1340     1448
##    m1[14]   3.21   3.21 0.24 0.23   2.82   3.60 1.00     1356     1595
##    m1[15]  -3.27  -3.27 0.32 0.31  -3.80  -2.76 1.00     1440     1547
##    m1[16]   0.46   0.46 0.17 0.16   0.17   0.74 1.00      762      948
##    m1[17]   2.79   2.79 0.22 0.21   2.44   3.16 1.00     1206     1460
##    m1[18]  -1.56  -1.57 0.23 0.22  -1.95  -1.19 1.00     1084     1097
##    m1[19]   3.32   3.32 0.25 0.23   2.92   3.72 1.00     1396     1595
##    m1[20]   3.66   3.66 0.26 0.25   3.24   4.09 1.00     1523     1631
##    m1[21]   2.20   2.20 0.20 0.19   1.88   2.53 1.00     1016     1395
##    m1[22]  -1.83  -1.82 0.30 0.30  -2.32  -1.35 1.00     1339     1448
##    m1[23]   1.57   1.58 0.18 0.17   1.28   1.88 1.00      859     1288
##    m1[24]   0.41   0.41 0.17 0.16   0.12   0.69 1.00      762      998
##    m1[25]   1.36   1.36 0.18 0.17   1.07   1.65 1.00      824     1271
##    m1[26]  -0.21  -0.21 0.24 0.24  -0.60   0.18 1.00     1063     1268
##    m1[27]   2.10   2.11 0.19 0.19   1.79   2.43 1.00      988     1343
##    m1[28]   0.97   0.96 0.22 0.22   0.61   1.34 1.00      975     1226
##    m1[29]  -2.21  -2.20 0.31 0.31  -2.73  -1.71 1.00     1412     1461
##    m1[30]   2.59   2.59 0.21 0.20   2.25   2.95 1.00     1139     1440
##    m1[31]   2.07   2.07 0.23 0.22   1.69   2.45 1.00     1029     1251
##    m1[32]   4.31   4.32 0.31 0.30   3.79   4.81 1.00     1413     1554
##    m1[33]   3.94   3.95 0.29 0.28   3.44   4.42 1.00     1337     1516
##    m1[34]   2.57   2.57 0.21 0.20   2.23   2.93 1.00     1131     1369
##    m1[35]   1.12   1.12 0.17 0.16   0.83   1.41 1.00      799     1115
##    m1[36]  -0.19  -0.19 0.18 0.17  -0.50   0.10 1.00      817     1025
##    m1[37]  -0.62  -0.62 0.19 0.18  -0.94  -0.30 1.00      893     1036
##    m1[38]  -1.51  -1.52 0.23 0.22  -1.90  -1.14 1.00     1073     1083
##    m1[39]   3.19   3.19 0.26 0.25   2.75   3.62 1.00     1194     1505
##    m1[40]  -3.00  -3.00 0.30 0.29  -3.50  -2.51 1.00     1386     1447
##    m1[41]   2.92   2.91 0.25 0.24   2.50   3.33 1.00     1146     1475
##    m1[42]   1.98   1.98 0.23 0.22   1.61   2.36 1.00     1020     1220
##    m1[43]   0.35   0.35 0.23 0.23  -0.02   0.73 1.00     1006     1233
##    m1[44]   3.93   3.93 0.28 0.27   3.49   4.39 1.00     1626     1624
##    m1[45]  -1.40  -1.40 0.22 0.21  -1.77  -1.03 1.00     1046      989
##    m1[46]   2.57   2.57 0.21 0.20   2.23   2.93 1.00     1131     1369
##    m1[47]  -0.50  -0.50 0.25 0.26  -0.91  -0.10 1.00     1101     1309
##    m1[48]   1.62   1.62 0.22 0.22   1.26   1.99 1.00      991     1238
##    m1[49]   2.89   2.89 0.23 0.21   2.53   3.26 1.00     1242     1482
##    m1[50]  -0.25  -0.25 0.18 0.17  -0.56   0.05 1.00      826     1003

fn(f1)
## Initial log joint probability = -99.5729 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -14.4515   8.86349e-06    0.00140858      0.4152      0.4152       29    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.31 -16.93 1.78 1.45 -20.65 -15.26 1.00      776      694
##    b[1]    -0.04  -0.04 0.19 0.19  -0.36   0.27 1.00     1243      991
##    b[2]     1.18   1.18 0.07 0.07   1.06   1.30 1.00     2515     1487
##    b[3]     0.45   0.44 0.27 0.25   0.03   0.91 1.00     1002      709
##    b[4]     0.87   0.87 0.24 0.24   0.47   1.27 1.00     1073      946
##    s        0.87   0.86 0.10 0.09   0.73   1.04 1.00     2343     1624
##    y1[1]    0.86   0.86 0.91 0.91  -0.58   2.31 1.00     1771     1817
##    y1[2]    0.63   0.64 0.89 0.89  -0.85   2.14 1.00     2054     1972
##    y1[3]    3.59   3.61 0.89 0.86   2.09   5.01 1.00     1662     1663
##    y1[4]    2.22   2.19 0.91 0.90   0.78   3.76 1.00     1757     1954
##    y1[5]   -2.31  -2.31 0.93 0.89  -3.85  -0.78 1.00     2015     1830
##    y1[6]   -1.89  -1.90 0.93 0.90  -3.40  -0.40 1.00     1982     1877
##    y1[7]    1.10   1.11 0.89 0.88  -0.32   2.59 1.00     1844     1811
##    y1[8]   -1.73  -1.75 0.88 0.83  -3.19  -0.24 1.00     1844     1877
##    y1[9]   -0.30  -0.30 0.92 0.90  -1.76   1.15 1.00     2100     1967
##    y1[10]   1.87   1.86 0.88 0.87   0.41   3.29 1.00     1851     1859
##    y1[11]   1.95   1.97 0.87 0.86   0.51   3.34 1.00     1977     1974
##    y1[12]   0.48   0.49 0.91 0.86  -1.01   1.97 1.00     2018     2014
##    y1[13]  -2.11  -2.11 0.91 0.90  -3.55  -0.60 1.00     2120     1882
##    y1[14]   2.72   2.73 0.91 0.88   1.16   4.20 1.00     1817     1895
##    y1[15]  -3.53  -3.54 0.91 0.91  -4.98  -2.04 1.00     1691     1808
##    y1[16]   0.06   0.07 0.92 0.90  -1.51   1.61 1.00     1837     1912
##    y1[17]   2.34   2.34 0.90 0.92   0.87   3.82 1.00     1603     1945
##    y1[18]  -1.00  -1.00 0.90 0.91  -2.43   0.47 1.00     1859     1943
##    y1[19]   2.83   2.82 0.91 0.91   1.36   4.31 1.00     2012     1959
##    y1[20]   4.04   4.01 0.91 0.90   2.57   5.57 1.00     1821     2046
##    y1[21]   1.77   1.75 0.90 0.87   0.34   3.31 1.00     2099     1689
##    y1[22]  -2.12  -2.10 0.91 0.88  -3.65  -0.67 1.00     1849     1852
##    y1[23]   2.04   2.03 0.90 0.87   0.51   3.46 1.00     1932     1939
##    y1[24]   0.01   0.02 0.88 0.89  -1.41   1.46 1.00     1931     1918
##    y1[25]   0.96   0.94 0.92 0.87  -0.55   2.44 1.00     1761     1820
##    y1[26]  -0.54  -0.54 0.92 0.89  -2.04   0.98 1.00     1797     1998
##    y1[27]   1.68   1.66 0.91 0.89   0.20   3.21 1.00     2076     1796
##    y1[28]   1.47   1.48 0.93 0.92  -0.06   3.00 1.00     1705     2012
##    y1[29]  -1.57  -1.55 0.93 0.92  -3.08  -0.06 1.00     2001     1955
##    y1[30]   3.03   3.02 0.88 0.86   1.60   4.49 1.00     2062     1845
##    y1[31]   2.51   2.51 0.92 0.92   1.00   4.03 1.00     1882     1969
##    y1[32]   3.84   3.81 0.93 0.92   2.28   5.38 1.00     1895     1833
##    y1[33]   4.35   4.37 0.91 0.91   2.86   5.81 1.00     1842     1933
##    y1[34]   2.98   2.97 0.89 0.85   1.54   4.43 1.00     2005     1970
##    y1[35]   1.59   1.57 0.89 0.89   0.19   3.07 1.00     2039     1932
##    y1[36]   0.33   0.33 0.92 0.90  -1.23   1.81 1.00     1886     1699
##    y1[37]  -0.10  -0.09 0.89 0.87  -1.51   1.39 1.00     1840     1974
##    y1[38]  -1.84  -1.82 0.90 0.89  -3.37  -0.38 1.00     2117     1871
##    y1[39]   3.65   3.66 0.91 0.91   2.16   5.17 1.00     1880     1915
##    y1[40]  -3.25  -3.24 0.91 0.86  -4.74  -1.78 1.00     2024     1897
##    y1[41]   3.37   3.36 0.92 0.92   1.90   4.90 1.00     1967     1931
##    y1[42]   2.41   2.41 0.90 0.90   0.93   3.89 1.00     2144     2020
##    y1[43]  -0.03  -0.03 0.90 0.88  -1.54   1.47 1.00     2038     2054
##    y1[44]   4.32   4.33 0.92 0.94   2.84   5.82 1.00     1853     1973
##    y1[45]  -1.74  -1.72 0.90 0.90  -3.22  -0.31 1.00     1941     2000
##    y1[46]   2.99   3.00 0.91 0.89   1.43   4.44 1.00     2021     1955
##    y1[47]  -0.83  -0.82 0.91 0.90  -2.31   0.68 1.00     2074     2100
##    y1[48]   1.25   1.24 0.92 0.87  -0.29   2.74 1.00     1807     1931
##    y1[49]   2.44   2.46 0.91 0.90   0.91   3.88 1.00     2033     1856
##    y1[50]  -0.57  -0.57 0.89 0.84  -2.06   0.86 1.00     1904     1745
##    m1[1]    0.87   0.87 0.20 0.19   0.54   1.19 1.00     1271      963
##    m1[2]    0.67   0.68 0.19 0.19   0.35   0.99 1.00     1254      994
##    m1[3]    3.59   3.59 0.27 0.26   3.16   4.04 1.00     1535     1312
##    m1[4]    2.26   2.26 0.26 0.26   1.85   2.69 1.00     1401     1175
##    m1[5]   -2.33  -2.32 0.32 0.31  -2.85  -1.81 1.00     1820     1623
##    m1[6]   -1.91  -1.91 0.28 0.29  -2.37  -1.45 1.00     1568     1468
##    m1[7]    1.10   1.10 0.21 0.21   0.75   1.45 1.00     1378     1238
##    m1[8]   -1.68  -1.68 0.22 0.21  -2.03  -1.31 1.00     1445     1087
##    m1[9]   -0.29  -0.28 0.24 0.23  -0.69   0.11 1.00     1524     1375
##    m1[10]   1.91   1.91 0.22 0.22   1.55   2.27 1.00     1429     1212
##    m1[11]   1.96   1.95 0.25 0.25   1.56   2.37 1.00     1361     1022
##    m1[12]   0.53   0.53 0.22 0.22   0.15   0.90 1.00     1436     1326
##    m1[13]  -2.11  -2.11 0.29 0.29  -2.57  -1.64 1.00     1603     1518
##    m1[14]   2.74   2.74 0.25 0.24   2.33   3.15 1.00     1595     1374
##    m1[15]  -3.52  -3.52 0.30 0.29  -3.99  -3.03 1.00     1806     1286
##    m1[16]   0.08   0.08 0.19 0.19  -0.23   0.39 1.00     1238      991
##    m1[17]   2.33   2.34 0.23 0.23   1.95   2.72 1.00     1511     1333
##    m1[18]  -1.00  -0.99 0.27 0.26  -1.42  -0.56 1.00     1620     1431
##    m1[19]   2.84   2.85 0.25 0.25   2.43   3.26 1.00     1616     1375
##    m1[20]   4.04   4.04 0.26 0.26   3.61   4.46 1.00     1650     1671
##    m1[21]   1.76   1.76 0.22 0.21   1.41   2.11 1.00     1400     1168
##    m1[22]  -2.10  -2.11 0.29 0.29  -2.56  -1.63 1.00     1602     1518
##    m1[23]   2.03   2.03 0.21 0.21   1.68   2.38 1.00     1364     1500
##    m1[24]   0.03   0.04 0.19 0.19  -0.28   0.35 1.00     1239      991
##    m1[25]   0.95   0.95 0.20 0.20   0.62   1.27 1.00     1279      973
##    m1[26]  -0.54  -0.55 0.25 0.25  -0.94  -0.14 1.00     1365     1096
##    m1[27]   1.67   1.67 0.21 0.21   1.33   2.02 1.00     1384     1087
##    m1[28]   1.47   1.47 0.25 0.24   1.07   1.86 1.00     1431     1126
##    m1[29]  -1.60  -1.60 0.33 0.30  -2.17  -1.07 1.00     1776     1602
##    m1[30]   3.02   3.01 0.23 0.23   2.64   3.39 1.00     1468     1578
##    m1[31]   2.53   2.53 0.25 0.24   2.13   2.93 1.00     1433     1242
##    m1[32]   3.83   3.83 0.31 0.31   3.32   4.35 1.00     1626     1517
##    m1[33]   4.34   4.33 0.29 0.28   3.87   4.83 1.00     1634     1313
##    m1[34]   2.99   2.99 0.23 0.23   2.62   3.36 1.00     1464     1584
##    m1[35]   1.59   1.59 0.21 0.21   1.24   1.93 1.00     1359     1413
##    m1[36]   0.32   0.33 0.23 0.22  -0.06   0.69 1.00     1454     1199
##    m1[37]  -0.08  -0.08 0.24 0.23  -0.48   0.30 1.00     1500     1380
##    m1[38]  -1.82  -1.82 0.23 0.22  -2.18  -1.44 1.00     1475     1167
##    m1[39]   3.62   3.61 0.27 0.26   3.18   4.06 1.00     1537     1312
##    m1[40]  -3.25  -3.25 0.29 0.28  -3.71  -2.78 1.00     1755     1251
##    m1[41]   3.35   3.35 0.26 0.25   2.93   3.78 1.00     1504     1260
##    m1[42]   2.45   2.45 0.25 0.24   2.05   2.85 1.00     1429     1252
##    m1[43]   0.00   0.00 0.24 0.24  -0.38   0.39 1.00     1315     1003
##    m1[44]   4.31   4.30 0.27 0.27   3.87   4.75 1.00     1693     1697
##    m1[45]  -1.71  -1.71 0.22 0.21  -2.06  -1.34 1.00     1450     1087
##    m1[46]   2.99   2.99 0.23 0.23   2.61   3.36 1.00     1464     1584
##    m1[47]  -0.82  -0.82 0.25 0.25  -1.23  -0.41 1.00     1396     1283
##    m1[48]   1.23   1.23 0.24 0.24   0.84   1.62 1.00     1305      937
##    m1[49]   2.43   2.44 0.24 0.24   2.04   2.82 1.00     1533     1351
##    m1[50]  -0.60  -0.60 0.20 0.19  -0.92  -0.28 1.00     1278      955

fn(f2)
## Initial log joint probability = -8920.61 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24       -12.484   9.52396e-06   0.000606877      0.2695      0.8802       51    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.91 -15.57 1.88 1.66 -19.37 -13.56 1.01      762      803
##    b[1]    -0.19  -0.18 0.20 0.20  -0.51   0.12 1.00      934     1094
##    b[2]     1.20   1.20 0.07 0.07   1.07   1.32 1.00     2107     1772
##    b[3]     0.85   0.84 0.33 0.33   0.34   1.38 1.00      691      718
##    b[4]     1.19   1.19 0.30 0.30   0.68   1.69 1.00      714      904
##    b[5]    -0.93  -0.92 0.51 0.48  -1.75  -0.09 1.00      482      550
##    s        0.84   0.84 0.09 0.09   0.71   1.01 1.00     2042     1703
##    y1[1]    0.73   0.73 0.87 0.85  -0.76   2.15 1.00     1842     1705
##    y1[2]    0.54   0.52 0.87 0.85  -0.88   2.00 1.00     1927     1900
##    y1[3]    3.28   3.30 0.90 0.89   1.79   4.74 1.00     1935     2013
##    y1[4]    2.57   2.54 0.89 0.92   1.08   4.04 1.00     1522     1773
##    y1[5]   -2.17  -2.16 0.90 0.89  -3.65  -0.73 1.00     2206     1906
##    y1[6]   -1.67  -1.68 0.89 0.86  -3.14  -0.18 1.00     2072     1883
##    y1[7]    1.28   1.31 0.90 0.89  -0.22   2.76 1.00     1785     1925
##    y1[8]   -1.88  -1.87 0.88 0.86  -3.35  -0.46 1.00     2028     1924
##    y1[9]   -0.10  -0.14 0.89 0.89  -1.53   1.35 1.00     1839     1816
##    y1[10]   1.82   1.82 0.87 0.88   0.39   3.23 1.00     1722     2016
##    y1[11]   2.23   2.22 0.90 0.89   0.70   3.70 1.00     1742     1851
##    y1[12]   0.68   0.65 0.89 0.87  -0.75   2.16 1.00     2095     1973
##    y1[13]  -1.90  -1.89 0.91 0.91  -3.38  -0.44 1.00     1728     1819
##    y1[14]   2.63   2.62 0.88 0.89   1.15   4.07 1.00     2116     2015
##    y1[15]  -3.72  -3.72 0.92 0.88  -5.19  -2.21 1.00     1775     1764
##    y1[16]  -0.06  -0.08 0.85 0.87  -1.43   1.37 1.00     1917     1846
##    y1[17]   2.22   2.22 0.86 0.83   0.79   3.61 1.00     1915     1885
##    y1[18]  -0.87  -0.86 0.88 0.85  -2.36   0.58 1.00     1888     1891
##    y1[19]   2.75   2.72 0.90 0.84   1.34   4.28 1.00     1752     1857
##    y1[20]   4.26   4.25 0.90 0.91   2.77   5.72 1.00     2008     1932
##    y1[21]   1.65   1.64 0.89 0.91   0.21   3.10 1.00     2061     1859
##    y1[22]  -1.90  -1.88 0.90 0.87  -3.37  -0.42 1.00     1629     1996
##    y1[23]   2.21   2.20 0.88 0.88   0.78   3.67 1.00     2111     1920
##    y1[24]  -0.09  -0.07 0.88 0.86  -1.50   1.38 1.00     1744     1895
##    y1[25]   0.81   0.81 0.87 0.84  -0.60   2.30 1.00     1836     1972
##    y1[26]  -0.30  -0.32 0.89 0.87  -1.73   1.15 1.00     1946     1968
##    y1[27]   1.54   1.56 0.88 0.86   0.11   2.96 1.00     1774     1940
##    y1[28]   1.12   1.13 0.91 0.88  -0.40   2.65 1.00     1810     1884
##    y1[29]  -2.02  -2.01 0.93 0.91  -3.56  -0.47 1.00     1721     1915
##    y1[30]   3.23   3.23 0.90 0.87   1.72   4.69 1.00     1666     1566
##    y1[31]   2.21   2.22 0.91 0.90   0.70   3.68 1.00     2120     1870
##    y1[32]   4.16   4.16 0.91 0.92   2.70   5.61 1.00     1550     1858
##    y1[33]   4.07   4.09 0.89 0.87   2.61   5.57 1.00     1905     1930
##    y1[34]   3.20   3.17 0.87 0.87   1.78   4.66 1.00     1905     1869
##    y1[35]   1.76   1.75 0.82 0.82   0.45   3.13 1.00     1990     1730
##    y1[36]   0.50   0.48 0.90 0.93  -0.93   2.00 1.00     1776     1766
##    y1[37]   0.09   0.08 0.85 0.83  -1.32   1.48 1.00     1743     1893
##    y1[38]  -2.00  -1.99 0.87 0.86  -3.44  -0.54 1.00     2016     2037
##    y1[39]   3.30   3.31 0.89 0.89   1.80   4.81 1.00     1936     2057
##    y1[40]  -3.45  -3.44 0.90 0.86  -4.92  -1.97 1.00     1728     2011
##    y1[41]   3.03   3.06 0.94 0.94   1.44   4.56 1.00     1773     1563
##    y1[42]   2.13   2.14 0.91 0.89   0.68   3.58 1.00     1705     1784
##    y1[43]   0.24   0.26 0.90 0.89  -1.25   1.69 1.00     1594     1772
##    y1[44]   4.56   4.54 0.90 0.91   3.07   6.03 1.00     1780     1930
##    y1[45]  -1.86  -1.86 0.88 0.87  -3.30  -0.38 1.00     1967     1873
##    y1[46]   3.20   3.22 0.87 0.88   1.75   4.59 1.00     1912     1899
##    y1[47]  -0.58  -0.56 0.88 0.87  -2.07   0.85 1.00     1881     1957
##    y1[48]   1.49   1.49 0.87 0.85   0.07   2.88 1.00     1892     1898
##    y1[49]   2.29   2.30 0.86 0.84   0.83   3.64 1.00     2028     1984
##    y1[50]  -0.77  -0.77 0.86 0.83  -2.18   0.57 1.00     1775     1808
##    m1[1]    0.74   0.74 0.20 0.20   0.41   1.06 1.00      972     1173
##    m1[2]    0.54   0.55 0.20 0.19   0.21   0.86 1.00      957     1165
##    m1[3]    3.28   3.29 0.31 0.30   2.79   3.77 1.00     1269     1422
##    m1[4]    2.55   2.54 0.30 0.29   2.08   3.06 1.00     1006     1094
##    m1[5]   -2.20  -2.20 0.32 0.33  -2.72  -1.68 1.00     1730     1604
##    m1[6]   -1.69  -1.69 0.30 0.29  -2.17  -1.19 1.00     1277     1509
##    m1[7]    1.29   1.29 0.23 0.23   0.90   1.66 1.00     1192     1322
##    m1[8]   -1.85  -1.84 0.23 0.23  -2.24  -1.47 1.00     1068     1431
##    m1[9]   -0.13  -0.12 0.25 0.25  -0.53   0.29 1.00     1368     1539
##    m1[10]   1.80   1.80 0.22 0.21   1.43   2.16 1.00     1109     1232
##    m1[11]   2.24   2.24 0.29 0.28   1.79   2.74 1.00      992     1115
##    m1[12]   0.70   0.70 0.24 0.23   0.31   1.10 1.00     1245     1371
##    m1[13]  -1.90  -1.89 0.30 0.30  -2.39  -1.38 1.00     1309     1488
##    m1[14]   2.64   2.64 0.25 0.25   2.22   3.06 1.00     1259     1407
##    m1[15]  -3.72  -3.71 0.31 0.30  -4.23  -3.20 1.00     1334     1767
##    m1[16]  -0.06  -0.06 0.20 0.19  -0.38   0.25 1.00      931     1132
##    m1[17]   2.23   2.23 0.23 0.23   1.83   2.62 1.00     1179     1248
##    m1[18]  -0.85  -0.84 0.27 0.28  -1.29  -0.40 1.00     1492     1623
##    m1[19]   2.75   2.75 0.25 0.25   2.32   3.17 1.00     1280     1521
##    m1[20]   4.28   4.27 0.29 0.28   3.79   4.76 1.00     1292     1516
##    m1[21]   1.65   1.65 0.22 0.21   1.28   2.00 1.00     1085     1208
##    m1[22]  -1.89  -1.88 0.30 0.30  -2.38  -1.38 1.00     1308     1488
##    m1[23]   2.23   2.23 0.24 0.24   1.84   2.62 1.00     1166     1312
##    m1[24]  -0.11  -0.10 0.20 0.19  -0.42   0.20 1.00      933     1106
##    m1[25]   0.82   0.83 0.20 0.20   0.50   1.15 1.00      979     1108
##    m1[26]  -0.30  -0.30 0.27 0.27  -0.74   0.15 1.00     1070     1167
##    m1[27]   1.56   1.56 0.21 0.20   1.19   1.91 1.00     1072     1206
##    m1[28]   1.12   1.12 0.31 0.29   0.61   1.63 1.00     1163     1198
##    m1[29]  -2.00  -1.99 0.40 0.38  -2.66  -1.34 1.00     1277     1255
##    m1[30]   3.23   3.23 0.26 0.25   2.82   3.66 1.00     1204     1358
##    m1[31]   2.20   2.21 0.30 0.29   1.71   2.70 1.00     1194     1228
##    m1[32]   4.14   4.14 0.35 0.35   3.58   4.73 1.00     1112     1287
##    m1[33]   4.04   4.06 0.32 0.32   3.52   4.55 1.00     1356     1463
##    m1[34]   3.21   3.21 0.26 0.25   2.79   3.63 1.00     1202     1360
##    m1[35]   1.78   1.78 0.23 0.23   1.40   2.16 1.00     1171     1270
##    m1[36]   0.50   0.50 0.24 0.23   0.10   0.89 1.00     1270     1401
##    m1[37]   0.08   0.08 0.25 0.25  -0.32   0.48 1.00     1334     1474
##    m1[38]  -1.99  -1.99 0.24 0.23  -2.39  -1.61 1.00     1088     1459
##    m1[39]   3.31   3.31 0.31 0.30   2.81   3.80 1.00     1271     1444
##    m1[40]  -3.45  -3.44 0.30 0.29  -3.94  -2.95 1.00     1302     1705
##    m1[41]   3.04   3.04 0.30 0.30   2.54   3.53 1.00     1248     1379
##    m1[42]   2.12   2.12 0.30 0.29   1.62   2.62 1.00     1191     1205
##    m1[43]   0.25   0.25 0.27 0.27  -0.18   0.69 1.00     1019     1072
##    m1[44]   4.54   4.54 0.30 0.28   4.05   5.04 1.00     1310     1539
##    m1[45]  -1.88  -1.88 0.23 0.23  -2.27  -1.50 1.00     1072     1387
##    m1[46]   3.21   3.20 0.26 0.25   2.79   3.63 1.00     1202     1360
##    m1[47]  -0.59  -0.59 0.27 0.27  -1.02  -0.13 1.00     1105     1226
##    m1[48]   1.50   1.50 0.28 0.27   1.06   1.96 1.00      978     1070
##    m1[49]   2.33   2.33 0.24 0.23   1.92   2.73 1.00     1198     1303
##    m1[50]  -0.75  -0.75 0.20 0.20  -1.09  -0.43 1.00      948     1185

tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,xa,y,col=ca),
             qplot(data=tb,xa,y,col=cb),
             qplot(data=tb,xb,y,col=ca),
             qplot(data=tb,xb,y,col=cb))

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)

fn(f0)
## Initial log joint probability = -127.979 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       28      -49.6612   0.000110023    0.00152751           1           1       33    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -52.38 -51.97 1.89 1.72 -55.90 -50.03 1.00      799     1344
##    b[1]    -0.05  -0.05 0.42 0.41  -0.74   0.65 1.00      761      832
##    b[2]     1.16   1.16 0.24 0.22   0.76   1.54 1.00     1996     1109
##    b[3]     1.18   1.18 0.22 0.22   0.82   1.54 1.00     2208     1446
##    b[4]     2.15   2.14 0.51 0.49   1.32   2.97 1.01      960      985
##    b[5]     1.79   1.80 0.51 0.49   0.93   2.63 1.00      745      982
##    s        1.78   1.76 0.20 0.18   1.49   2.14 1.00     1916     1374
##    y1[1]    5.17   5.12 1.85 1.82   2.16   8.29 1.00     1778     1893
##    y1[2]    3.92   3.91 1.86 1.83   0.95   7.00 1.00     1844     1802
##    y1[3]    4.77   4.77 1.93 1.96   1.59   8.02 1.00     1824     1860
##    y1[4]    3.76   3.74 1.91 1.93   0.59   6.85 1.00     1929     1974
##    y1[5]    4.58   4.56 1.87 1.79   1.44   7.59 1.00     1988     1977
##    y1[6]    0.09   0.07 1.85 1.84  -2.99   3.13 1.00     1947     1941
##    y1[7]    5.98   5.99 1.88 1.87   2.89   9.14 1.00     1895     2012
##    y1[8]    2.09   2.12 1.84 1.76  -0.95   5.14 1.00     1828     1980
##    y1[9]   -1.08  -1.06 1.85 1.85  -4.16   1.93 1.00     1830     2012
##    y1[10]   2.30   2.34 1.84 1.79  -0.82   5.35 1.00     2023     1981
##    y1[11]   4.77   4.78 1.93 1.91   1.51   7.91 1.00     1938     1936
##    y1[12]   2.09   2.06 1.90 1.90  -1.00   5.24 1.00     2075     1848
##    y1[13]   4.10   4.11 1.90 1.88   0.96   7.10 1.00     1776     1893
##    y1[14]   4.33   4.31 1.94 1.90   1.24   7.52 1.00     1793     1742
##    y1[15]   1.15   1.11 1.84 1.79  -1.77   4.21 1.00     1928     1939
##    y1[16]   0.95   0.95 1.89 1.87  -2.17   4.04 1.00     2018     2054
##    y1[17]   1.68   1.72 1.85 1.79  -1.46   4.68 1.00     1739     1718
##    y1[18]   1.16   1.13 1.89 1.84  -1.88   4.31 1.00     1862     1842
##    y1[19]  -0.68  -0.68 1.90 1.85  -3.76   2.40 1.00     1957     1895
##    y1[20]   4.15   4.15 1.81 1.79   1.22   7.07 1.00     1875     1746
##    y1[21]   0.54   0.52 1.87 1.83  -2.52   3.57 1.00     2049     1872
##    y1[22]   3.01   3.01 1.86 1.82   0.10   6.12 1.00     1709     1704
##    y1[23]  -0.48  -0.42 1.81 1.76  -3.62   2.41 1.00     1816     1917
##    y1[24]   3.00   2.97 1.83 1.85   0.03   5.82 1.00     2058     1811
##    y1[25]   4.21   4.22 1.85 1.86   1.24   7.23 1.00     2073     1918
##    y1[26]   2.00   1.99 1.87 1.88  -1.00   5.14 1.00     1918     1864
##    y1[27]   1.06   1.05 1.90 1.86  -2.05   4.13 1.00     1860     1666
##    y1[28]   0.17   0.18 1.86 1.93  -2.90   3.17 1.00     1739     1841
##    y1[29]  -1.10  -1.06 1.87 1.83  -4.12   1.91 1.00     1838     1946
##    y1[30]  -1.92  -1.92 1.92 1.87  -5.01   1.21 1.00     1751     1790
##    y1[31]   2.30   2.35 1.87 1.86  -0.85   5.43 1.00     1956     1864
##    y1[32]  -0.94  -0.93 1.86 1.80  -4.04   2.10 1.00     2010     1895
##    y1[33]  -1.87  -1.87 1.89 1.92  -4.94   1.19 1.00     1954     1744
##    y1[34]  -0.61  -0.60 1.86 1.81  -3.65   2.48 1.00     1788     1665
##    y1[35]   2.26   2.24 1.89 1.99  -0.80   5.30 1.00     2035     1917
##    y1[36]  -0.73  -0.73 1.91 1.88  -3.87   2.42 1.00     1895     1974
##    y1[37]  -2.15  -2.13 1.95 1.88  -5.41   1.05 1.00     1890     1857
##    y1[38]   3.60   3.58 1.83 1.78   0.67   6.67 1.00     2116     2044
##    y1[39]   1.29   1.34 1.84 1.84  -1.81   4.34 1.00     1785     1893
##    y1[40]   0.09   0.05 1.85 1.89  -2.84   3.12 1.00     1940     1971
##    y1[41]   0.35   0.30 1.90 1.88  -2.76   3.48 1.00     2054     1973
##    y1[42]   3.65   3.60 1.83 1.75   0.60   6.70 1.00     1956     2015
##    y1[43]   3.64   3.65 1.88 1.83   0.46   6.61 1.00     1955     1805
##    y1[44]   0.78   0.80 1.81 1.78  -2.21   3.77 1.00     1900     2014
##    y1[45]   1.05   1.08 1.89 1.84  -2.10   4.12 1.00     1898     1856
##    y1[46]   2.02   2.05 1.84 1.79  -0.99   5.10 1.00     1855     1892
##    y1[47]  -1.24  -1.20 1.91 1.95  -4.28   1.79 1.00     1706     1823
##    y1[48]   1.30   1.30 1.85 1.86  -1.76   4.42 1.00     1813     1967
##    y1[49]   3.95   3.96 1.89 1.87   0.78   7.04 1.00     1846     1796
##    y1[50]   1.25   1.28 1.86 1.84  -1.71   4.19 1.00     1644     2014
##    m1[1]    5.18   5.18 0.55 0.53   4.26   6.07 1.01     1460     1528
##    m1[2]    3.86   3.85 0.55 0.52   2.93   4.74 1.00     1598     1548
##    m1[3]    4.79   4.77 0.61 0.60   3.78   5.79 1.00     1473     1471
##    m1[4]    3.66   3.65 0.61 0.59   2.64   4.67 1.00     1537     1406
##    m1[5]    4.51   4.51 0.49 0.48   3.69   5.30 1.00     1407     1726
##    m1[6]    0.14   0.14 0.53 0.53  -0.74   1.02 1.00     1822     1524
##    m1[7]    5.93   5.93 0.64 0.61   4.85   6.96 1.00     1541     1362
##    m1[8]    2.14   2.14 0.52 0.51   1.28   3.01 1.00     1034     1166
##    m1[9]   -1.07  -1.07 0.59 0.59  -2.08  -0.11 1.00     1907     1719
##    m1[10]   2.27   2.26 0.49 0.49   1.51   3.09 1.00      984     1244
##    m1[11]   4.80   4.80 0.62 0.60   3.79   5.80 1.00     2135     1539
##    m1[12]   2.06   2.08 0.72 0.69   0.90   3.28 1.00     1934     1542
##    m1[13]   4.14   4.15 0.67 0.64   3.00   5.21 1.00     1528     1430
##    m1[14]   4.30   4.33 0.69 0.66   3.19   5.41 1.00     1707     1463
##    m1[15]   1.22   1.21 0.50 0.48   0.39   2.06 1.00      954     1178
##    m1[16]   0.92   0.91 0.55 0.54   0.02   1.84 1.00     1953     1362
##    m1[17]   1.67   1.68 0.54 0.53   0.77   2.57 1.00     1369     1417
##    m1[18]   1.14   1.15 0.65 0.62   0.07   2.21 1.00     1341     1340
##    m1[19]  -0.67  -0.67 0.59 0.58  -1.64   0.27 1.00     1610     1551
##    m1[20]   4.14   4.15 0.54 0.52   3.23   4.99 1.00     1409     1313
##    m1[21]   0.55   0.55 0.61 0.61  -0.45   1.56 1.00     1402     1389
##    m1[22]   3.00   2.98 0.57 0.57   2.09   3.93 1.00     1158     1321
##    m1[23]  -0.48  -0.48 0.54 0.52  -1.39   0.39 1.00     1476     1355
##    m1[24]   3.01   3.00 0.48 0.46   2.23   3.80 1.00     1275     1412
##    m1[25]   4.26   4.25 0.48 0.47   3.47   5.03 1.00     1345     1311
##    m1[26]   1.99   1.97 0.51 0.49   1.17   2.85 1.00     1654     1290
##    m1[27]   1.05   1.06 0.56 0.55   0.13   1.98 1.00     1818     1634
##    m1[28]   0.18   0.18 0.51 0.51  -0.68   1.02 1.00     1526     1361
##    m1[29]  -1.11  -1.12 0.62 0.62  -2.12  -0.08 1.00     1199     1362
##    m1[30]  -1.93  -1.92 0.69 0.65  -3.08  -0.84 1.00     1858     1629
##    m1[31]   2.26   2.26 0.57 0.56   1.32   3.21 1.00     1406     1385
##    m1[32]  -0.98  -0.98 0.59 0.58  -1.99  -0.04 1.00     1890     1688
##    m1[33]  -1.84  -1.84 0.55 0.53  -2.74  -0.94 1.00     1035     1342
##    m1[34]  -0.67  -0.67 0.61 0.60  -1.67   0.32 1.00     1653     1511
##    m1[35]   2.27   2.26 0.49 0.49   1.50   3.08 1.00      981     1181
##    m1[36]  -0.68  -0.68 0.68 0.68  -1.80   0.42 1.00     1144     1396
##    m1[37]  -2.12  -2.11 0.69 0.68  -3.28  -1.01 1.00     1584     1304
##    m1[38]   3.65   3.65 0.50 0.49   2.85   4.48 1.00     1784     1694
##    m1[39]   1.35   1.34 0.48 0.47   0.57   2.14 1.00      909     1166
##    m1[40]   0.10   0.10 0.50 0.49  -0.75   0.94 1.00     1589     1571
##    m1[41]   0.32   0.32 0.66 0.65  -0.74   1.45 1.00     1316     1257
##    m1[42]   3.64   3.63 0.46 0.44   2.89   4.39 1.00     1369     1187
##    m1[43]   3.67   3.66 0.53 0.51   2.81   4.56 1.00     1398     1385
##    m1[44]   0.79   0.80 0.47 0.46   0.01   1.58 1.00     1439     1431
##    m1[45]   1.05   1.05 0.58 0.56   0.10   2.01 1.00     1055     1183
##    m1[46]   2.06   2.06 0.45 0.44   1.32   2.79 1.00     1557     1566
##    m1[47]  -1.27  -1.26 0.61 0.59  -2.27  -0.28 1.00     1515     1385
##    m1[48]   1.33   1.33 0.43 0.44   0.64   2.06 1.00      830     1031
##    m1[49]   3.95   3.94 0.55 0.53   3.06   4.87 1.00     1626     1520
##    m1[50]   1.33   1.32 0.48 0.47   0.55   2.12 1.00      908     1180

fn(f1)
## Initial log joint probability = -125.882 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       35      -49.6599   0.000110695     0.0017254      0.4031           1       41    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -53.07 -52.70 2.13 1.94 -56.99 -50.29 1.00      544      945
##    b[1]    -0.03  -0.03 0.54 0.54  -0.92   0.88 1.00      687      880
##    b[2]     1.16   1.16 0.24 0.24   0.76   1.56 1.00     2298     1450
##    b[3]     1.17   1.17 0.24 0.24   0.78   1.56 1.00     1659     1536
##    b[4]     2.07   2.07 0.82 0.80   0.74   3.40 1.01      581      778
##    b[5]     1.77   1.78 0.76 0.75   0.51   3.01 1.01      640      847
##    b[6]     0.07   0.02 1.17 1.06  -1.83   1.93 1.01      548      765
##    s        1.80   1.79 0.20 0.20   1.51   2.16 1.00     2078     1346
##    y1[1]    5.09   5.11 1.94 1.97   1.92   8.23 1.00     1975     1806
##    y1[2]    3.84   3.87 1.93 1.85   0.61   6.87 1.00     1934     1969
##    y1[3]    4.75   4.75 1.97 1.98   1.42   7.97 1.00     1734     1700
##    y1[4]    3.71   3.67 1.96 1.98   0.53   6.82 1.00     1962     1916
##    y1[5]    4.50   4.48 1.88 1.86   1.50   7.57 1.00     2105     1840
##    y1[6]    0.07   0.06 1.90 1.85  -3.15   3.22 1.00     1804     1909
##    y1[7]    5.97   6.00 1.92 1.92   2.77   9.15 1.00     2062     2015
##    y1[8]    2.23   2.28 1.85 1.83  -0.83   5.23 1.00     2169     1990
##    y1[9]   -1.09  -1.10 1.91 1.94  -4.23   1.95 1.00     1893     1794
##    y1[10]   2.32   2.31 1.88 1.87  -0.76   5.43 1.00     1961     2009
##    y1[11]   4.82   4.77 2.03 2.06   1.58   8.18 1.00     2007     1973
##    y1[12]   2.03   2.02 1.95 1.85  -1.29   5.22 1.00     2172     1915
##    y1[13]   4.07   4.07 1.99 1.83   0.86   7.43 1.00     1998     2088
##    y1[14]   4.21   4.16 1.92 1.88   1.10   7.38 1.00     2106     1894
##    y1[15]   1.24   1.20 1.91 1.86  -1.88   4.35 1.00     1657     1804
##    y1[16]   0.87   0.85 1.95 1.95  -2.33   4.12 1.00     1875     1885
##    y1[17]   1.68   1.71 1.92 1.96  -1.58   4.71 1.00     1849     2000
##    y1[18]   1.09   1.14 1.92 1.87  -2.05   4.26 1.00     1977     1855
##    y1[19]  -0.69  -0.65 1.91 1.95  -3.78   2.43 1.00     1972     1807
##    y1[20]   4.12   4.08 1.91 1.89   0.97   7.28 1.00     2057     1916
##    y1[21]   0.56   0.60 1.90 1.80  -2.63   3.63 1.00     1755     1909
##    y1[22]   3.04   3.05 1.89 1.83  -0.05   6.19 1.00     1966     1919
##    y1[23]  -0.54  -0.54 1.89 1.88  -3.64   2.43 1.00     2144     1836
##    y1[24]   2.92   2.92 1.93 1.93  -0.25   6.01 1.00     1722     1820
##    y1[25]   4.25   4.29 1.89 1.83   1.08   7.35 1.00     1985     1820
##    y1[26]   1.94   2.00 1.95 1.86  -1.33   5.14 1.00     1919     1838
##    y1[27]   1.03   1.03 1.93 1.97  -2.14   4.15 1.00     1741     1917
##    y1[28]   0.15   0.17 1.87 1.89  -2.96   3.16 1.00     1770     1788
##    y1[29]  -1.04  -0.99 1.90 1.81  -4.18   1.99 1.00     1714     1725
##    y1[30]  -1.96  -1.92 1.97 1.96  -5.12   1.26 1.00     2026     1815
##    y1[31]   2.30   2.18 1.97 1.94  -0.82   5.74 1.00     2011     1922
##    y1[32]  -1.03  -0.92 1.96 1.94  -4.25   2.06 1.00     1910     1836
##    y1[33]  -1.84  -1.85 1.94 1.92  -4.96   1.49 1.00     1500     1682
##    y1[34]  -0.73  -0.73 1.89 1.80  -3.98   2.34 1.00     2073     1949
##    y1[35]   2.25   2.20 1.87 1.85  -0.76   5.30 1.00     1809     1762
##    y1[36]  -0.67  -0.69 1.98 1.98  -3.89   2.63 1.00     2035     1969
##    y1[37]  -2.11  -2.07 1.93 1.91  -5.34   0.93 1.00     2126     1901
##    y1[38]   3.56   3.55 1.92 1.86   0.38   6.67 1.00     1912     1894
##    y1[39]   1.31   1.35 1.89 1.89  -1.75   4.37 1.00     1890     1782
##    y1[40]   0.09   0.12 1.86 1.81  -2.97   3.11 1.00     1785     1985
##    y1[41]   0.37   0.41 1.97 1.93  -2.95   3.60 1.00     1560     1861
##    y1[42]   3.59   3.59 1.85 1.81   0.50   6.58 1.00     1994     2057
##    y1[43]   3.62   3.62 1.93 1.92   0.42   6.72 1.00     1660     1840
##    y1[44]   0.74   0.75 1.88 1.89  -2.41   3.74 1.00     2134     2009
##    y1[45]   1.10   1.10 1.98 1.89  -2.22   4.28 1.00     1960     2016
##    y1[46]   2.07   2.07 1.88 1.89  -0.97   5.21 1.00     1886     1875
##    y1[47]  -1.22  -1.20 1.91 1.91  -4.36   1.96 1.00     2050     1910
##    y1[48]   1.38   1.35 1.92 1.88  -1.75   4.57 1.00     1983     1809
##    y1[49]   3.93   3.94 1.90 1.87   0.82   7.09 1.00     1879     1800
##    y1[50]   1.36   1.37 1.87 1.83  -1.74   4.44 1.00     1643     1854
##    m1[1]    5.16   5.16 0.62 0.61   4.16   6.19 1.00     2063     1256
##    m1[2]    3.84   3.85 0.64 0.61   2.80   4.86 1.00     2045     1196
##    m1[3]    4.73   4.72 0.76 0.73   3.46   5.97 1.00      979     1167
##    m1[4]    3.67   3.67 0.69 0.66   2.53   4.78 1.00     1549     1480
##    m1[5]    4.50   4.49 0.58 0.56   3.55   5.44 1.00     1904     1660
##    m1[6]    0.10   0.11 0.55 0.55  -0.79   1.01 1.00     1582     1482
##    m1[7]    5.91   5.91 0.70 0.68   4.79   7.05 1.00     2200     1220
##    m1[8]    2.15   2.14 0.57 0.55   1.22   3.12 1.00     1296     1542
##    m1[9]   -1.11  -1.12 0.60 0.60  -2.12  -0.13 1.00     1863     1664
##    m1[10]   2.29   2.28 0.54 0.53   1.37   3.18 1.00     1164     1303
##    m1[11]   4.78   4.78 0.78 0.79   3.52   6.01 1.00     1137     1354
##    m1[12]   2.03   2.02 0.74 0.71   0.81   3.25 1.00     1715     1375
##    m1[13]   4.14   4.12 0.77 0.76   2.87   5.41 1.00     1729     1857
##    m1[14]   4.27   4.27 0.76 0.73   3.03   5.52 1.00     2317     1207
##    m1[15]   1.23   1.22 0.56 0.55   0.27   2.19 1.00     1057     1384
##    m1[16]   0.92   0.92 0.73 0.69  -0.30   2.12 1.00     1244     1415
##    m1[17]   1.68   1.67 0.58 0.57   0.73   2.63 1.00     1820     1615
##    m1[18]   1.07   1.06 0.74 0.72  -0.12   2.27 1.00     1490     1327
##    m1[19]  -0.68  -0.67 0.61 0.60  -1.77   0.31 1.00     1868     1625
##    m1[20]   4.13   4.12 0.65 0.65   3.08   5.19 1.00     1580     1723
##    m1[21]   0.57   0.57 0.62 0.61  -0.45   1.59 1.00     2019     1681
##    m1[22]   3.02   3.02 0.61 0.58   1.98   4.04 1.00     1520     1376
##    m1[23]  -0.48  -0.46 0.55 0.54  -1.43   0.42 1.00     1914     1382
##    m1[24]   2.96   2.97 0.60 0.57   1.97   3.91 1.00      950     1075
##    m1[25]   4.25   4.24 0.58 0.59   3.32   5.16 1.00     1632     1634
##    m1[26]   1.99   1.99 0.67 0.64   0.86   3.07 1.00     1272     1494
##    m1[27]   1.01   1.03 0.59 0.58   0.08   1.96 1.00     1528     1146
##    m1[28]   0.18   0.18 0.55 0.54  -0.79   1.06 1.00     1632     1630
##    m1[29]  -1.08  -1.09 0.78 0.74  -2.37   0.22 1.00      842     1067
##    m1[30]  -1.98  -1.98 0.70 0.70  -3.11  -0.83 1.00     2096     1943
##    m1[31]   2.26   2.26 0.61 0.61   1.25   3.29 1.00     1819     1467
##    m1[32]  -1.03  -1.03 0.60 0.60  -2.02  -0.06 1.00     1836     1664
##    m1[33]  -1.81  -1.82 0.72 0.70  -2.99  -0.63 1.00      696     1009
##    m1[34]  -0.68  -0.66 0.64 0.62  -1.81   0.35 1.00     1855     1666
##    m1[35]   2.28   2.28 0.54 0.53   1.38   3.18 1.00     1177     1397
##    m1[36]  -0.68  -0.68 0.75 0.71  -1.96   0.61 1.00     1114     1382
##    m1[37]  -2.11  -2.10 0.69 0.67  -3.28  -1.01 1.00     2044     1572
##    m1[38]   3.64   3.66 0.63 0.63   2.61   4.63 1.00     1139     1320
##    m1[39]   1.36   1.35 0.54 0.53   0.44   2.28 1.00     1003     1423
##    m1[40]   0.05   0.06 0.54 0.52  -0.85   0.91 1.00     1464     1437
##    m1[41]   0.36   0.37 0.77 0.74  -0.94   1.64 1.00     1056     1313
##    m1[42]   3.63   3.62 0.58 0.58   2.68   4.55 1.00     1480     1641
##    m1[43]   3.62   3.63 0.65 0.60   2.52   4.71 1.00      976     1104
##    m1[44]   0.74   0.76 0.52 0.49  -0.13   1.57 1.00     1250     1397
##    m1[45]   1.06   1.04 0.63 0.62   0.01   2.13 1.00     1244     1447
##    m1[46]   2.05   2.06 0.54 0.54   1.15   2.90 1.00     1235     1411
##    m1[47]  -1.26  -1.24 0.61 0.58  -2.30  -0.29 1.00     2048     1600
##    m1[48]   1.35   1.34 0.51 0.51   0.49   2.21 1.00      862     1044
##    m1[49]   3.95   3.97 0.66 0.64   2.86   5.03 1.00     1279     1462
##    m1[50]   1.34   1.33 0.54 0.53   0.42   2.26 1.00      999     1423

fn(f2)
## Initial log joint probability = -101.671 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -12.1861   6.24919e-05   0.000665711       0.796       0.796       32    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.76 -16.46 2.09 1.97 -20.59 -14.01 1.00      831     1308
##    b[1]    -0.16  -0.17 0.25 0.25  -0.57   0.23 1.00      623      840
##    b[2]     0.68   0.68 0.12 0.12   0.47   0.87 1.00     2306     1607
##    b[3]     1.28   1.29 0.11 0.11   1.10   1.46 1.00     1312     1545
##    b[4]     2.58   2.60 0.38 0.37   1.92   3.19 1.01      513      743
##    b[5]     2.39   2.40 0.36 0.35   1.77   2.96 1.01      592      717
##    b[6]    -1.24  -1.24 0.10 0.11  -1.41  -1.07 1.00     2235     1707
##    b[7]    -1.39  -1.41 0.53 0.53  -2.26  -0.49 1.01      477      695
##    s        0.86   0.85 0.09 0.09   0.72   1.02 1.00     1744     1386
##    y1[1]    4.82   4.83 0.89 0.87   3.36   6.27 1.00     1937     1931
##    y1[2]    4.84   4.80 0.92 0.90   3.37   6.34 1.00     1950     1919
##    y1[3]    3.01   3.01 0.97 0.96   1.46   4.60 1.00     1725     1958
##    y1[4]    3.62   3.61 0.91 0.87   2.10   5.12 1.00     1995     1766
##    y1[5]    4.10   4.09 0.88 0.87   2.66   5.54 1.00     1896     2058
##    y1[6]    0.29   0.29 0.90 0.85  -1.17   1.79 1.00     1699     1657
##    y1[7]    5.55   5.56 0.91 0.92   4.10   7.06 1.00     2084     1946
##    y1[8]    1.32   1.33 0.92 0.89  -0.26   2.81 1.00     2168     1940
##    y1[9]   -2.54  -2.52 0.90 0.90  -4.01  -1.11 1.00     2111     1733
##    y1[10]   0.49   0.50 0.91 0.89  -0.99   2.00 1.00     1735     1793
##    y1[11]   3.63   3.64 0.95 0.91   2.07   5.14 1.00     1749     1806
##    y1[12]   5.72   5.70 0.98 0.97   4.14   7.35 1.00     1829     1874
##    y1[13]   5.50   5.50 0.95 0.94   3.96   7.03 1.00     1992     1934
##    y1[14]   7.29   7.29 0.97 0.96   5.70   8.83 1.00     1906     1707
##    y1[15]   1.53   1.53 0.89 0.88   0.10   2.98 1.00     1843     1973
##    y1[16]  -1.17  -1.16 0.93 0.94  -2.69   0.34 1.00     1867     1947
##    y1[17]   2.73   2.73 0.91 0.95   1.24   4.19 1.00     1678     1875
##    y1[18]   3.70   3.69 0.98 0.97   2.14   5.31 1.00     1621     1545
##    y1[19]  -0.37  -0.38 0.90 0.91  -1.80   1.07 1.00     2079     1930
##    y1[20]   4.05   4.04 0.91 0.89   2.56   5.54 1.00     1959     1896
##    y1[21]   1.93   1.91 0.92 0.92   0.47   3.49 1.00     1957     1718
##    y1[22]   0.32   0.32 0.94 0.95  -1.20   1.85 1.00     1841     2017
##    y1[23]  -0.67  -0.69 0.90 0.88  -2.19   0.81 1.00     1914     1746
##    y1[24]   2.95   2.95 0.92 0.90   1.45   4.45 1.00     1920     1931
##    y1[25]   3.62   3.61 0.91 0.90   2.18   5.15 1.00     1804     1997
##    y1[26]   1.11   1.10 0.94 0.90  -0.45   2.64 1.00     1889     1870
##    y1[27]   2.38   2.38 0.88 0.89   0.93   3.82 1.00     2101     2039
##    y1[28]   0.91   0.91 0.92 0.93  -0.66   2.41 1.00     2095     1704
##    y1[29]   0.39   0.37 0.96 0.90  -1.19   2.01 1.00     1610     1675
##    y1[30]  -4.53  -4.54 0.95 0.93  -6.09  -2.92 1.00     1853     1791
##    y1[31]   3.38   3.39 0.92 0.91   1.78   4.87 1.00     1800     1856
##    y1[32]  -2.28  -2.27 0.90 0.86  -3.78  -0.81 1.00     1971     1859
##    y1[33]  -2.12  -2.13 0.91 0.89  -3.67  -0.64 1.00     1848     1780
##    y1[34]   0.07   0.10 0.92 0.91  -1.47   1.55 1.00     1621     1702
##    y1[35]   0.56   0.57 0.91 0.92  -0.94   2.08 1.00     1858     1845
##    y1[36]   2.51   2.51 0.98 0.97   0.88   4.04 1.00     2032     2015
##    y1[37]  -4.48  -4.46 0.97 0.98  -6.11  -2.91 1.00     1858     1737
##    y1[38]   3.29   3.33 0.92 0.87   1.77   4.76 1.00     1830     1949
##    y1[39]   1.19   1.19 0.87 0.90  -0.21   2.60 1.00     1552     1817
##    y1[40]  -0.22  -0.22 0.88 0.86  -1.67   1.23 1.00     1838     1844
##    y1[41]   2.57   2.59 0.94 0.91   0.99   4.10 1.00     1608     1879
##    y1[42]   3.11   3.07 0.90 0.90   1.67   4.59 1.00     1871     1820
##    y1[43]   3.21   3.17 0.92 0.91   1.72   4.74 1.00     1896     2059
##    y1[44]   0.87   0.89 0.87 0.84  -0.56   2.33 1.00     1954     1876
##    y1[45]   2.63   2.64 0.89 0.88   1.20   4.13 1.00     1803     1972
##    y1[46]   2.90   2.92 0.92 0.90   1.37   4.37 1.00     2061     1875
##    y1[47]  -2.40  -2.39 0.90 0.90  -3.86  -0.91 1.00     2064     2037
##    y1[48]   0.50   0.53 0.90 0.88  -1.00   1.95 1.00     1669     1821
##    y1[49]   2.90   2.90 0.91 0.92   1.45   4.40 1.00     1753     1858
##    y1[50]   1.19   1.20 0.91 0.92  -0.31   2.65 1.00     1705     1808
##    m1[1]    4.81   4.81 0.28 0.29   4.36   5.28 1.00     2095     1727
##    m1[2]    4.82   4.82 0.30 0.30   4.33   5.31 1.00     2252     1571
##    m1[3]    3.03   3.04 0.39 0.39   2.39   3.67 1.00     1258     1371
##    m1[4]    3.59   3.59 0.32 0.32   3.08   4.13 1.00     1254     1281
##    m1[5]    4.10   4.10 0.26 0.27   3.66   4.53 1.00     2022     1553
##    m1[6]    0.29   0.29 0.26 0.25  -0.14   0.72 1.00     1719     1833
##    m1[7]    5.55   5.55 0.31 0.32   5.02   6.04 1.00     2180     1502
##    m1[8]    1.33   1.32 0.27 0.27   0.89   1.77 1.00     1106     1505
##    m1[9]   -2.53  -2.52 0.31 0.31  -3.04  -2.00 1.00     2145     1787
##    m1[10]   0.52   0.51 0.29 0.29   0.06   1.00 1.00     1125     1444
##    m1[11]   3.64   3.64 0.37 0.38   3.02   4.25 1.00     1242     1370
##    m1[12]   5.73   5.73 0.48 0.47   4.98   6.49 1.00     1694     1540
##    m1[13]   5.52   5.52 0.38 0.38   4.90   6.15 1.00     1933     1589
##    m1[14]   7.26   7.26 0.42 0.42   6.57   7.93 1.00     2524     1499
##    m1[15]   1.56   1.56 0.26 0.27   1.14   1.98 1.00      953     1209
##    m1[16]  -1.13  -1.13 0.38 0.39  -1.78  -0.50 1.00     1212     1632
##    m1[17]   2.74   2.74 0.29 0.29   2.27   3.22 1.00     1227     1368
##    m1[18]   3.69   3.69 0.40 0.39   3.03   4.34 1.00     1139     1321
##    m1[19]  -0.35  -0.35 0.30 0.30  -0.84   0.15 1.00     1810     1562
##    m1[20]   4.04   4.04 0.29 0.30   3.54   4.52 1.00     1700     1557
##    m1[21]   1.93   1.94 0.32 0.31   1.41   2.45 1.00     1383     1494
##    m1[22]   0.33   0.32 0.36 0.37  -0.24   0.92 1.00     1396     1385
##    m1[23]  -0.66  -0.66 0.28 0.27  -1.12  -0.20 1.00     1646     1649
##    m1[24]   2.92   2.93 0.28 0.27   2.46   3.39 1.00      949     1191
##    m1[25]   3.63   3.63 0.27 0.27   3.18   4.06 1.00     1749     1506
##    m1[26]   1.10   1.10 0.32 0.31   0.57   1.61 1.00     1292     1530
##    m1[27]   2.38   2.38 0.30 0.28   1.89   2.87 1.00     1514     1566
##    m1[28]   0.90   0.91 0.28 0.27   0.46   1.34 1.00     1426     1527
##    m1[29]   0.38   0.38 0.38 0.35  -0.24   1.01 1.00      879     1100
##    m1[30]  -4.55  -4.54 0.41 0.41  -5.22  -3.85 1.00     2252     1911
##    m1[31]   3.36   3.36 0.31 0.31   2.86   3.87 1.00     1224     1368
##    m1[32]  -2.33  -2.32 0.31 0.30  -2.83  -1.81 1.00     2117     1823
##    m1[33]  -2.13  -2.13 0.33 0.33  -2.65  -1.58 1.00      620      995
##    m1[34]   0.06   0.07 0.32 0.32  -0.47   0.59 1.00     1818     1650
##    m1[35]   0.58   0.58 0.28 0.28   0.13   1.06 1.00     1120     1401
##    m1[36]   2.51   2.50 0.43 0.43   1.80   3.21 1.00     1816     1561
##    m1[37]  -4.46  -4.46 0.41 0.39  -5.13  -3.78 1.00     2024     1493
##    m1[38]   3.26   3.25 0.30 0.30   2.77   3.74 1.00     1085     1345
##    m1[39]   1.18   1.17 0.25 0.26   0.78   1.59 1.00      866     1401
##    m1[40]  -0.23  -0.22 0.26 0.25  -0.66   0.18 1.00     1614     1591
##    m1[41]   2.58   2.57 0.41 0.40   1.89   3.26 1.00     1268     1462
##    m1[42]   3.12   3.13 0.27 0.26   2.68   3.55 1.00     1608     1611
##    m1[43]   3.22   3.23 0.31 0.30   2.71   3.72 1.00     1052     1304
##    m1[44]   0.87   0.88 0.25 0.24   0.47   1.27 1.00     1294     1520
##    m1[45]   2.62   2.62 0.31 0.33   2.11   3.15 1.00     1502     1379
##    m1[46]   2.89   2.89 0.26 0.26   2.46   3.31 1.00     1009     1121
##    m1[47]  -2.39  -2.39 0.32 0.31  -2.93  -1.85 1.00     1902     1572
##    m1[48]   0.52   0.51 0.24 0.25   0.13   0.92 1.00      779     1186
##    m1[49]   2.95   2.95 0.32 0.32   2.44   3.48 1.00     1288     1441
##    m1[50]   1.19   1.18 0.25 0.26   0.79   1.59 1.00      863     1401

generalized linear regression

log normal regression

objective variable [0,Infinity)

# of samples is N,  
log mi=sum(bj*xji),j=[0,K],i=[1,N]  
log yi~N(mi,s)  

ex5-2.stan

log normal regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~lognormal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,log(x1+x2),1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex5-2.stan')

mle=mdl$optimize(data=data)  # it sometimes occur error and stop process
## Initial log joint probability = -3273.07 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       36      -19.7607   8.31624e-05    0.00157544           1           1       57    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -19.76
##    b[1]      -1.22
##    b[2]       0.23
##    b[3]       0.15
##    s          0.65
##    y1[1]      0.92
##    y1[2]      2.79
##    y1[3]      0.55
##    y1[4]      2.40
##    y1[5]      2.62
##    y1[6]      0.90
##    y1[7]      1.96
##    y1[8]      0.11
##    y1[9]      3.75
##    y1[10]     0.95
##    y1[11]     1.92
##    y1[12]     2.02
##    y1[13]     5.38
##    y1[14]     1.31
##    y1[15]     0.80
##    y1[16]     0.92
##    y1[17]     0.65
##    y1[18]     2.13
##    y1[19]     1.81
##    y1[20]     2.48
##    m1[1]     -0.74
##    m1[2]      1.00
##    m1[3]     -0.12
##    m1[4]      1.20
##    m1[5]      0.68
##    m1[6]      0.21
##    m1[7]      0.56
##    m1[8]     -0.76
##    m1[9]      1.94
##    m1[10]     0.58
##    m1[11]     0.33
##    m1[12]     0.40
##    m1[13]     1.04
##    m1[14]     1.23
##    m1[15]     1.28
##    m1[16]     0.62
##    m1[17]    -0.07
##    m1[18]     1.22
##    m1[19]     0.49
##    m1[20]     0.20
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.34 -21.98  1.51 1.36 -25.37 -20.57 1.00      555      846
##    b[1]    -1.24  -1.24  0.49 0.46  -2.04  -0.43 1.00      312      371
##    b[2]     0.23   0.23  0.07 0.07   0.11   0.35 1.00      469      630
##    b[3]     0.16   0.16  0.07 0.06   0.04   0.28 1.00      971      841
##    s        0.76   0.74  0.14 0.13   0.57   1.01 1.00      861     1077
##    y1[1]    0.72   0.47  0.81 0.36   0.11   2.11 1.00     1074     1196
##    y1[2]    3.96   2.72  4.49 2.06   0.73  11.14 1.00     1936     1821
##    y1[3]    1.25   0.88  1.29 0.63   0.23   3.32 1.00     1362     1697
##    y1[4]    4.68   3.35  5.29 2.41   0.88  12.81 1.00     1867     1778
##    y1[5]    2.68   1.99  2.52 1.42   0.53   7.25 1.00     2324     2098
##    y1[6]    1.77   1.22  2.10 0.89   0.31   4.81 1.00     1727     1918
##    y1[7]    2.42   1.76  2.42 1.24   0.50   6.25 1.00     1924     1834
##    y1[8]    0.65   0.45  0.72 0.33   0.11   1.69 1.00     1032     1660
##    y1[9]   10.86   7.48 12.96 5.68   1.68  31.40 1.00     1677     1772
##    y1[10]   2.41   1.70  2.64 1.22   0.49   6.27 1.00     1912     1599
##    y1[11]   1.98   1.42  2.05 1.10   0.36   5.43 1.00     1719     1799
##    y1[12]   2.09   1.54  2.19 1.08   0.42   5.48 1.00     1664     1843
##    y1[13]   4.09   2.98  4.28 2.09   0.76  11.04 1.00     1960     1931
##    y1[14]   5.22   3.61  5.69 2.73   0.89  14.19 1.00     1887     2057
##    y1[15]   5.10   3.63  5.38 2.63   0.96  14.00 1.00     2221     1890
##    y1[16]   2.67   1.95  2.87 1.45   0.50   7.00 1.00     2047     1968
##    y1[17]   1.37   0.92  1.60 0.73   0.24   3.80 1.00     1519     1682
##    y1[18]   4.76   3.44  4.60 2.60   0.86  12.84 1.00     1741     1785
##    y1[19]   2.47   1.64  3.19 1.25   0.39   7.07 1.00     2039     1859
##    y1[20]   1.72   1.24  1.99 0.90   0.32   4.67 1.00     1933     1882
##    m1[1]   -0.75  -0.76  0.38 0.35  -1.40  -0.12 1.00      321      382
##    m1[2]    1.01   1.01  0.27 0.26   0.59   1.45 1.00     1942     1462
##    m1[3]   -0.12  -0.12  0.25 0.23  -0.54   0.29 1.00      374      492
##    m1[4]    1.21   1.22  0.28 0.27   0.75   1.67 1.00     1342     1382
##    m1[5]    0.68   0.69  0.26 0.26   0.24   1.11 1.00     2059     1694
##    m1[6]    0.21   0.22  0.29 0.28  -0.26   0.69 1.00     1088      959
##    m1[7]    0.57   0.57  0.17 0.17   0.29   0.86 1.00     1396     1215
##    m1[8]   -0.77  -0.77  0.39 0.36  -1.42  -0.13 1.00      321      391
##    m1[9]    1.96   1.96  0.42 0.41   1.28   2.64 1.00      730      958
##    m1[10]   0.59   0.59  0.17 0.17   0.31   0.87 1.00     1534     1208
##    m1[11]   0.33   0.33  0.28 0.27  -0.13   0.78 1.00     1377     1193
##    m1[12]   0.41   0.41  0.18 0.17   0.12   0.70 1.00      803      822
##    m1[13]   1.05   1.05  0.21 0.22   0.71   1.40 1.00     1938     1461
##    m1[14]   1.25   1.24  0.30 0.30   0.77   1.75 1.00     1757     1493
##    m1[15]   1.29   1.29  0.25 0.25   0.88   1.71 1.00     1264     1539
##    m1[16]   0.63   0.63  0.32 0.31   0.09   1.13 1.00     1983     1333
##    m1[17]  -0.06  -0.07  0.37 0.35  -0.68   0.52 1.00      687      791
##    m1[18]   1.23   1.23  0.30 0.29   0.72   1.73 1.00     1367     1475
##    m1[19]   0.50   0.50  0.35 0.34  -0.08   1.08 1.00     1521     1197
##    m1[20]   0.20   0.20  0.23 0.22  -0.18   0.58 1.00      689      805

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

poisson regression

objective variable [0,Infinity), integer. variance of error is near to mean  
(normal linear regression, correlation is near to 1,-1, variance of error is constant)  

# of samples is N,  
log li=sum(bj*xji),j=[0,K],i=[1,N]  
yi~Po(li)  
 or  
li=sum(bj*xji),j=[0,k]  
yi~Po(exp li)  

when xj -> xj +1, y -> y* exp bj   

ex6-1.stan

poisson regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] l=X*b;
  y~poisson_log(l);
}
generated quantities{
  array[N] int y1;
  vector[N] l1=X*b;
  for(i in 1:N){
    y1[i]=poisson_log_rng(l1[i]);
  }
}
n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
          y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='poisson')
## 
## Call:  glm(formula = f, family = "poisson", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##      -0.223        1.635        0.547  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       59.9 
## Residual Deviance: 22    AIC: 81.5
X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -242.607 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13      -4.98764   5.28446e-05   0.000456263      0.8447      0.8447       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -4.99
##    b[1]      -0.22
##    b[2]       1.64
##    b[3]       0.55
##    y1[1]      1.00
##    y1[2]      0.00
##    y1[3]      0.00
##    y1[4]      1.00
##    y1[5]      0.00
##    y1[6]      3.00
##    y1[7]      0.00
##    y1[8]      0.00
##    y1[9]      1.00
##    y1[10]     1.00
##    y1[11]     4.00
##    y1[12]     2.00
##    y1[13]     8.00
##    y1[14]     0.00
##    y1[15]     2.00
##    y1[16]     2.00
##    y1[17]     1.00
##    y1[18]     0.00
##    y1[19]     3.00
##    y1[20]     0.00
##    y1[21]     5.00
##    y1[22]     8.00
##    y1[23]     0.00
##    y1[24]     0.00
##    y1[25]     1.00
##    y1[26]     1.00
##    y1[27]     0.00
##    y1[28]     0.00
##    y1[29]     4.00
##    y1[30]     1.00
##    l1[1]      0.42
##    l1[2]      0.42
##    l1[3]     -0.97
##    l1[4]      0.10
##    l1[5]      0.21
##    l1[6]      0.62
##    l1[7]      0.35
##    l1[8]      0.38
##    l1[9]     -0.93
##    l1[10]     0.72
##    l1[11]     1.12
##    l1[12]     1.16
##    l1[13]     1.90
##    l1[14]    -1.08
##    l1[15]     0.38
##    l1[16]    -0.64
##    l1[17]    -0.43
##    l1[18]    -0.93
##    l1[19]     1.33
##    l1[20]    -0.75
##    l1[21]     1.84
##    l1[22]     1.39
##    l1[23]    -1.12
##    l1[24]     0.21
##    l1[25]    -0.70
##    l1[26]    -0.33
##    l1[27]    -0.55
##    l1[28]    -0.33
##    l1[29]     0.88
##    l1[30]    -0.32
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='l1',ch=F)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   -6.44  -6.14 1.19 0.96 -8.80 -5.15 1.00      721     1058
##    b[1]   -0.26  -0.25 0.27 0.27 -0.72  0.14 1.00      916      905
##    b[2]    1.65   1.64 0.29 0.30  1.18  2.14 1.00     1300     1192
##    b[3]    0.55   0.54 0.29 0.28  0.08  1.03 1.00      888      839
##    y1[1]   1.52   1.00 1.26 1.48  0.00  4.00 1.00     2041     1804
##    y1[2]   1.54   1.00 1.28 1.48  0.00  4.00 1.00     1608     1747
##    y1[3]   0.38   0.00 0.62 0.00  0.00  2.00 1.00     1975     1828
##    y1[4]   1.10   1.00 1.12 1.48  0.00  3.00 1.00     2011     1763
##    y1[5]   1.19   1.00 1.15 1.48  0.00  3.00 1.00     1831     1824
##    y1[6]   1.81   2.00 1.37 1.48  0.00  4.00 1.00     1930     1896
##    y1[7]   1.45   1.00 1.24 1.48  0.00  4.00 1.00     1677     1903
##    y1[8]   1.45   1.00 1.22 1.48  0.00  4.00 1.00     2029     2014
##    y1[9]   0.40   0.00 0.65 0.00  0.00  2.00 1.00     1897     1906
##    y1[10]  2.03   2.00 1.47 1.48  0.00  5.00 1.00     2148     2108
##    y1[11]  3.07   3.00 1.91 1.48  0.00  6.00 1.00     1960     1989
##    y1[12]  3.17   3.00 1.84 1.48  1.00  6.00 1.00     1772     1858
##    y1[13]  6.60   6.00 3.00 2.97  2.00 12.00 1.00     1955     1832
##    y1[14]  0.34   0.00 0.59 0.00  0.00  1.00 1.00     1769     1797
##    y1[15]  1.48   1.00 1.25 1.48  0.00  4.00 1.01     1866     1867
##    y1[16]  0.52   0.00 0.74 0.00  0.00  2.00 1.00     1877     1938
##    y1[17]  0.66   0.00 0.84 0.00  0.00  2.00 1.00     1860     1901
##    y1[18]  0.42   0.00 0.67 0.00  0.00  2.00 1.00     1954     2006
##    y1[19]  3.76   3.00 2.22 1.48  1.00  8.00 1.00     1874     1655
##    y1[20]  0.47   0.00 0.70 0.00  0.00  2.00 1.00     2015     1920
##    y1[21]  6.35   6.00 2.79 2.97  2.00 11.00 1.00     1780     1997
##    y1[22]  3.93   4.00 2.18 1.48  1.00  8.00 1.00     2083     1962
##    y1[23]  0.33   0.00 0.60 0.00  0.00  1.00 1.00     1948     2002
##    y1[24]  1.22   1.00 1.12 1.48  0.00  3.00 1.00     1819     1720
##    y1[25]  0.52   0.00 0.77 0.00  0.00  2.00 1.00     2107     1797
##    y1[26]  0.72   0.00 0.89 0.00  0.00  2.00 1.00     1792     1714
##    y1[27]  0.58   0.00 0.82 0.00  0.00  2.00 1.00     1991     1883
##    y1[28]  0.70   1.00 0.84 1.48  0.00  2.00 1.00     1805     1648
##    y1[29]  2.40   2.00 1.63 1.48  0.00  5.00 1.00     1770     1922
##    y1[30]  0.74   1.00 0.87 1.48  0.00  2.00 1.00     1939     1882
##    l1[1]   0.38   0.38 0.21 0.22  0.02  0.72 1.00     1591     1521
##    l1[2]   0.39   0.39 0.23 0.23 -0.02  0.76 1.00      989     1077
##    l1[3]  -1.02  -1.00 0.35 0.36 -1.60 -0.48 1.00      949      851
##    l1[4]   0.06   0.08 0.24 0.24 -0.35  0.44 1.00      935      984
##    l1[5]   0.17   0.18 0.24 0.24 -0.23  0.54 1.00      951      948
##    l1[6]   0.58   0.58 0.20 0.20  0.26  0.89 1.00     1698     1635
##    l1[7]   0.32   0.33 0.23 0.23 -0.08  0.69 1.00      976     1058
##    l1[8]   0.35   0.35 0.22 0.22 -0.02  0.69 1.00     1573     1546
##    l1[9]  -0.97  -0.96 0.35 0.36 -1.55 -0.44 1.00      946      820
##    l1[10]  0.69   0.69 0.19 0.19  0.38  0.98 1.00     1758     1596
##    l1[11]  1.09   1.09 0.26 0.25  0.64  1.50 1.00     1183     1266
##    l1[12]  1.13   1.13 0.17 0.17  0.84  1.40 1.00     1862     1473
##    l1[13]  1.87   1.88 0.22 0.22  1.52  2.22 1.00     1542     1619
##    l1[14] -1.12  -1.10 0.37 0.38 -1.73 -0.56 1.00      959      874
##    l1[15]  0.35   0.35 0.22 0.22 -0.03  0.69 1.00     1573     1546
##    l1[16] -0.68  -0.67 0.31 0.32 -1.20 -0.20 1.00      922      806
##    l1[17] -0.47  -0.46 0.29 0.29 -0.96 -0.03 1.00      908      798
##    l1[18] -0.98  -0.97 0.41 0.43 -1.65 -0.33 1.01     1345     1202
##    l1[19]  1.30   1.31 0.28 0.27  0.82  1.74 1.01     1204     1380
##    l1[20] -0.79  -0.78 0.32 0.33 -1.34 -0.30 1.00      930      801
##    l1[21]  1.81   1.82 0.21 0.21  1.46  2.15 1.00     1577     1628
##    l1[22]  1.36   1.36 0.18 0.18  1.06  1.64 1.00     1816     1593
##    l1[23] -1.17  -1.15 0.44 0.46 -1.89 -0.47 1.01     1337     1249
##    l1[24]  0.18   0.18 0.24 0.24 -0.23  0.56 1.00     1510     1472
##    l1[25] -0.74  -0.73 0.37 0.38 -1.35 -0.15 1.01     1362     1197
##    l1[26] -0.37  -0.36 0.31 0.32 -0.90  0.14 1.01     1398     1367
##    l1[27] -0.59  -0.58 0.35 0.36 -1.17 -0.03 1.01     1373     1289
##    l1[28] -0.37  -0.37 0.28 0.28 -0.84  0.05 1.00      906      905
##    l1[29]  0.84   0.85 0.24 0.24  0.42  1.23 1.00     1128     1179
##    l1[30] -0.36  -0.35 0.31 0.32 -0.88  0.14 1.01     1400     1367

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4 6
##   0 7 3 0 0 0 0
##   1 4 3 0 0 0 0
##   2 0 3 1 0 0 0
##   3 0 1 1 1 0 1
##   4 0 0 1 1 1 0
##   5 0 0 0 1 0 0
##   7 0 0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4 6
##   0 3 3 0 0 0 0
##   1 3 0 0 0 0 0
##   2 0 1 0 0 0 0
##   3 0 1 1 1 0 0
##   4 0 0 0 1 0 0
##   5 0 0 0 0 0 0
##   7 0 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4 6
##   0 4 0 0 0 0 0
##   1 1 3 0 0 0 0
##   2 0 2 1 0 0 0
##   3 0 0 0 0 0 1
##   4 0 0 1 0 1 0
##   5 0 0 0 1 0 0
##   7 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 6
##   0 9 1 0 0 0
##   1 5 2 0 0 0
##   2 0 4 0 0 0
##   3 0 1 1 1 1
##   4 0 1 0 2 0
##   5 0 0 0 1 0
##   7 0 0 0 0 1
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 6
##   0 5 1 0 0 0
##   1 3 0 0 0 0
##   2 0 1 0 0 0
##   3 0 1 1 1 0
##   4 0 0 0 1 0
##   5 0 0 0 0 0
##   7 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 6
##   0 4 0 0 0 0
##   1 2 2 0 0 0
##   2 0 3 0 0 0
##   3 0 0 0 0 1
##   4 0 1 0 1 0
##   5 0 0 0 1 0
##   7 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

logistic regression

# of samples is N,  
objective variable 0/1 binary  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~Ber(pi), 0/1 binary  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-2.stan

logistic regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~bernoulli_logit(z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=bernoulli_rng(inv_logit(z1[i]));
  }
}
n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='binomial') # it can caluculte when all trials are once
## 
## Call:  glm(formula = f, family = "binomial", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.201        0.167        0.324  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       41.1 
## Residual Deviance: 40.9  AIC: 46.9
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -38.9076 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -20.4263   0.000312115   7.48949e-05      0.8677      0.8677       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -20.43
##    b[1]       0.20
##    b[2]       0.17
##    b[3]       0.32
##    y1[1]      1.00
##    y1[2]      0.00
##    y1[3]      1.00
##    y1[4]      1.00
##    y1[5]      0.00
##    y1[6]      0.00
##    y1[7]      1.00
##    y1[8]      1.00
##    y1[9]      1.00
##    y1[10]     1.00
##    y1[11]     1.00
##    y1[12]     0.00
##    y1[13]     1.00
##    y1[14]     0.00
##    y1[15]     0.00
##    y1[16]     1.00
##    y1[17]     1.00
##    y1[18]     1.00
##    y1[19]     0.00
##    y1[20]     1.00
##    y1[21]     1.00
##    y1[22]     1.00
##    y1[23]     1.00
##    y1[24]     1.00
##    y1[25]     1.00
##    y1[26]     1.00
##    y1[27]     0.00
##    y1[28]     1.00
##    y1[29]     1.00
##    y1[30]     0.00
##    z1[1]      0.35
##    z1[2]      0.18
##    z1[3]      0.53
##    z1[4]      0.17
##    z1[5]      0.10
##    z1[6]      0.27
##    z1[7]      0.16
##    z1[8]      0.09
##    z1[9]      0.31
##    z1[10]     0.21
##    z1[11]     0.07
##    z1[12]     0.23
##    z1[13]     0.59
##    z1[14]     0.62
##    z1[15]     0.12
##    z1[16]     0.52
##    z1[17]     0.28
##    z1[18]     0.12
##    z1[19]     0.19
##    z1[20]     0.20
##    z1[21]     0.42
##    z1[22]     0.56
##    z1[23]     0.15
##    z1[24]     0.40
##    z1[25]     0.08
##    z1[26]     0.15
##    z1[27]     0.13
##    z1[28]     0.11
##    z1[29]     0.37
##    z1[30]     0.46
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.00 -21.68 1.26 1.03 -24.55 -20.61 1.00      917     1112
##    b[1]     0.22   0.23 0.47 0.46  -0.54   1.02 1.00     1639     1365
##    b[2]     0.17   0.13 0.81 0.78  -1.13   1.56 1.00     1388     1220
##    b[3]     0.36   0.36 0.91 0.91  -1.07   1.88 1.00     1083      970
##    y1[1]    0.58   1.00 0.49 0.00   0.00   1.00 1.00     1846       NA
##    y1[2]    0.53   1.00 0.50 0.00   0.00   1.00 1.00     2016       NA
##    y1[3]    0.64   1.00 0.48 0.00   0.00   1.00 1.00     1982       NA
##    y1[4]    0.54   1.00 0.50 0.00   0.00   1.00 1.00     1859       NA
##    y1[5]    0.52   1.00 0.50 0.00   0.00   1.00 1.00     1963       NA
##    y1[6]    0.57   1.00 0.50 0.00   0.00   1.00 1.00     2047       NA
##    y1[7]    0.55   1.00 0.50 0.00   0.00   1.00 1.00     1954       NA
##    y1[8]    0.52   1.00 0.50 0.00   0.00   1.00 1.00     1958       NA
##    y1[9]    0.56   1.00 0.50 0.00   0.00   1.00 1.00     1820       NA
##    y1[10]   0.56   1.00 0.50 0.00   0.00   1.00 1.00     1874       NA
##    y1[11]   0.51   1.00 0.50 0.00   0.00   1.00 1.00     1939       NA
##    y1[12]   0.57   1.00 0.50 0.00   0.00   1.00 1.00     1964       NA
##    y1[13]   0.64   1.00 0.48 0.00   0.00   1.00 1.00     1811       NA
##    y1[14]   0.65   1.00 0.48 0.00   0.00   1.00 1.00     1733       NA
##    y1[15]   0.53   1.00 0.50 0.00   0.00   1.00 1.00     1855       NA
##    y1[16]   0.62   1.00 0.48 0.00   0.00   1.00 1.00     1833       NA
##    y1[17]   0.55   1.00 0.50 0.00   0.00   1.00 1.00     1841       NA
##    y1[18]   0.53   1.00 0.50 0.00   0.00   1.00 1.00     1973       NA
##    y1[19]   0.55   1.00 0.50 0.00   0.00   1.00 1.00     1932       NA
##    y1[20]   0.56   1.00 0.50 0.00   0.00   1.00 1.00     1894       NA
##    y1[21]   0.60   1.00 0.49 0.00   0.00   1.00 1.00     1879       NA
##    y1[22]   0.63   1.00 0.48 0.00   0.00   1.00 1.00     1746       NA
##    y1[23]   0.53   1.00 0.50 0.00   0.00   1.00 1.00     1934       NA
##    y1[24]   0.61   1.00 0.49 0.00   0.00   1.00 1.00     2047       NA
##    y1[25]   0.52   1.00 0.50 0.00   0.00   1.00 1.00     2039       NA
##    y1[26]   0.52   1.00 0.50 0.00   0.00   1.00 1.00     1605       NA
##    y1[27]   0.56   1.00 0.50 0.00   0.00   1.00 1.00     2089       NA
##    y1[28]   0.53   1.00 0.50 0.00   0.00   1.00 1.00     2022       NA
##    y1[29]   0.59   1.00 0.49 0.00   0.00   1.00 1.00     1847       NA
##    y1[30]   0.62   1.00 0.48 0.00   0.00   1.00 1.00     1954       NA
##    z1[1]    0.37   0.35 0.95 0.94  -1.12   1.97 1.00     1389     1219
##    z1[2]    0.20   0.20 0.45 0.43  -0.55   0.95 1.00     1711     1389
##    z1[3]    0.59   0.57 0.77 0.75  -0.64   1.90 1.00     1401     1542
##    z1[4]    0.18   0.19 0.45 0.43  -0.57   0.93 1.00     1731     1330
##    z1[5]    0.11   0.13 0.58 0.58  -0.81   1.08 1.00     1670     1518
##    z1[6]    0.28   0.28 0.63 0.61  -0.70   1.31 1.00     1471     1328
##    z1[7]    0.17   0.18 0.46 0.44  -0.58   0.94 1.00     1747     1366
##    z1[8]    0.10   0.12 0.62 0.62  -0.87   1.12 1.00     1636     1518
##    z1[9]    0.33   0.30 0.79 0.77  -0.90   1.64 1.00     1424     1340
##    z1[10]   0.23   0.25 0.49 0.48  -0.57   1.06 1.00     1596     1519
##    z1[11]   0.08   0.10 0.68 0.69  -1.00   1.21 1.00     1580     1536
##    z1[12]   0.25   0.26 0.52 0.52  -0.60   1.13 1.00     1550     1406
##    z1[13]   0.65   0.63 0.85 0.82  -0.72   2.05 1.00     1403     1370
##    z1[14]   0.68   0.68 0.92 0.92  -0.78   2.14 1.00     1399     1364
##    z1[15]   0.13   0.14 0.53 0.52  -0.71   1.00 1.00     1709     1471
##    z1[16]   0.58   0.56 0.77 0.75  -0.63   1.89 1.00     1402     1335
##    z1[17]   0.30   0.28 0.68 0.66  -0.77   1.42 1.00     1449     1320
##    z1[18]   0.13   0.14 0.53 0.52  -0.71   1.00 1.00     1709     1471
##    z1[19]   0.21   0.22 0.46 0.44  -0.55   0.98 1.00     1672     1430
##    z1[20]   0.22   0.23 0.47 0.46  -0.56   1.03 1.00     1628     1449
##    z1[21]   0.47   0.43 0.89 0.89  -0.98   1.95 1.00     1410     1365
##    z1[22]   0.62   0.59 0.81 0.79  -0.65   1.95 1.00     1401     1501
##    z1[23]   0.17   0.17 0.47 0.45  -0.59   0.93 1.00     1751     1406
##    z1[24]   0.45   0.42 0.95 0.95  -1.08   2.01 1.00     1409     1322
##    z1[25]   0.09   0.11 0.66 0.66  -0.95   1.18 1.00     1600     1540
##    z1[26]   0.16   0.17 0.47 0.45  -0.60   0.94 1.00     1753     1406
##    z1[27]   0.14   0.14 0.51 0.49  -0.68   0.97 1.00     1723     1519
##    z1[28]   0.13   0.13 0.55 0.54  -0.76   1.04 1.00     1692     1448
##    z1[29]   0.39   0.35 1.02 1.01  -1.24   2.11 1.00     1383     1203
##    z1[30]   0.52   0.48 0.80 0.79  -0.78   1.85 1.00     1402     1405

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##      1
##   0 13
##   1 17
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##      1
##   0 10
##   1 12
## 
## , ,  = b
## 
##    
##      1
##   0  3
##   1  5
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##      1
##   0 13
##   1 17
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##      1
##   0 10
##   1 12
## 
## , ,  = b
## 
##    map
##      1
##   0  3
##   1  5
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

multi logistic regression

# of samples is N,  
# of trials of a sample i is mi,  
objective variable [0,n], integer  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~B(mi, pi), # of acutual incident  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-3.stan

multi logistic regression

data{
  int N;
  int K;
  array[N] int m;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~binomial_logit(m,z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
  }
}
n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)

mdl=cmdstan_model('./ex6-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -155.274 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -96.8953   0.000390871    8.9511e-05           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -96.90
##    b[1]       0.07
##    b[2]       0.88
##    b[3]       0.99
##    y1[1]      2.00
##    y1[2]      2.00
##    y1[3]      0.00
##    y1[4]      3.00
##    y1[5]      3.00
##    y1[6]      0.00
##    y1[7]      7.00
##    y1[8]      4.00
##    y1[9]      5.00
##    y1[10]     5.00
##    y1[11]     1.00
##    y1[12]     6.00
##    y1[13]     7.00
##    y1[14]     2.00
##    y1[15]     2.00
##    y1[16]     4.00
##    y1[17]     0.00
##    y1[18]     2.00
##    y1[19]     1.00
##    y1[20]     2.00
##    y1[21]     4.00
##    y1[22]     3.00
##    y1[23]     6.00
##    y1[24]     4.00
##    y1[25]     3.00
##    y1[26]     2.00
##    y1[27]     3.00
##    y1[28]     3.00
##    y1[29]     4.00
##    y1[30]     2.00
##    z1[1]      1.48
##    z1[2]      0.84
##    z1[3]     -0.34
##    z1[4]      1.83
##    z1[5]     -0.34
##    z1[6]     -0.50
##    z1[7]      1.38
##    z1[8]      0.51
##    z1[9]      1.29
##    z1[10]    -0.58
##    z1[11]    -0.40
##    z1[12]     0.65
##    z1[13]     1.87
##    z1[14]    -0.67
##    z1[15]    -0.16
##    z1[16]     0.42
##    z1[17]    -0.71
##    z1[18]     0.31
##    z1[19]     0.29
##    z1[20]     0.92
##    z1[21]     0.29
##    z1[22]     0.51
##    z1[23]     1.81
##    z1[24]     0.90
##    z1[25]     1.22
##    z1[26]    -0.73
##    z1[27]    -0.67
##    z1[28]     1.03
##    z1[29]     0.58
##    z1[30]     0.53
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -98.46 -98.11 1.32 1.03 -101.16 -97.07 1.00      694     1140
##    b[1]     0.08   0.07 0.33 0.33   -0.45   0.65 1.00      624      925
##    b[2]     0.91   0.91 0.37 0.36    0.32   1.54 1.00      908     1081
##    b[3]     1.01   1.01 0.37 0.37    0.42   1.62 1.00      682      943
##    y1[1]    2.43   3.00 0.69 0.00    1.00   3.00 1.00     1994       NA
##    y1[2]    1.39   1.00 0.66 1.48    0.00   2.00 1.00     2110       NA
##    y1[3]    0.41   0.00 0.49 0.00    0.00   1.00 1.00     1849       NA
##    y1[4]    2.57   3.00 0.62 0.00    1.00   3.00 1.00     1787       NA
##    y1[5]    2.51   2.00 1.25 1.48    1.00   5.00 1.00     1893     1772
##    y1[6]    0.74   1.00 0.69 1.48    0.00   2.00 1.00     1912       NA
##    y1[7]    7.19   7.00 1.27 1.48    5.00   9.00 1.00     1844       NA
##    y1[8]    4.36   4.00 1.40 1.48    2.00   7.00 1.00     1520       NA
##    y1[9]    4.73   5.00 1.02 1.48    3.00   6.00 1.00     1924       NA
##    y1[10]   3.22   3.00 1.57 1.48    1.00   6.00 1.00     1646     1746
##    y1[11]   1.17   1.00 0.87 1.48    0.00   3.00 1.00     2021       NA
##    y1[12]   4.64   5.00 1.28 1.48    3.00   7.00 1.00     1970       NA
##    y1[13]   6.89   7.00 1.07 1.48    5.00   8.00 1.00     1882       NA
##    y1[14]   3.06   3.00 1.52 1.48    1.00   6.00 1.00     1823     1814
##    y1[15]   4.18   4.00 1.67 1.48    2.00   7.00 1.00     1479     1832
##    y1[16]   5.49   6.00 1.58 1.48    3.00   8.00 1.00     1942     2088
##    y1[17]   0.31   0.00 0.46 0.00    0.00   1.00 1.00     1962       NA
##    y1[18]   1.14   1.00 0.71 1.48    0.00   2.00 1.00     1893       NA
##    y1[19]   1.72   2.00 0.87 1.48    0.00   3.00 1.00     1953       NA
##    y1[20]   1.43   2.00 0.65 0.00    0.00   2.00 1.00     2025       NA
##    y1[21]   2.80   3.00 1.18 1.48    1.00   5.00 1.00     1453       NA
##    y1[22]   3.75   4.00 1.23 1.48    2.00   6.00 1.00     1765       NA
##    y1[23]   5.17   5.00 0.88 1.48    4.00   6.00 1.00     1708       NA
##    y1[24]   4.27   4.00 1.11 1.48    2.00   6.00 1.00     1816       NA
##    y1[25]   2.30   2.00 0.75 1.48    1.00   3.00 1.00     1876       NA
##    y1[26]   2.00   2.00 1.22 1.48    0.00   4.00 1.00     1931     1848
##    y1[27]   1.98   2.00 1.21 1.48    0.00   4.00 1.00     1947     1932
##    y1[28]   5.20   5.00 1.17 1.48    3.00   7.00 1.00     1846       NA
##    y1[29]   5.16   5.00 1.44 1.48    3.00   7.00 1.00     2117     1799
##    y1[30]   3.80   4.00 1.23 1.48    2.00   6.00 1.00     2144       NA
##    z1[1]    1.52   1.52 0.36 0.34    0.95   2.10 1.00     1238     1102
##    z1[2]    0.86   0.85 0.23 0.22    0.50   1.25 1.00     1704     1540
##    z1[3]   -0.35  -0.35 0.28 0.27   -0.79   0.12 1.00      780      821
##    z1[4]    1.89   1.88 0.48 0.46    1.11   2.68 1.00     1110     1151
##    z1[5]   -0.34  -0.35 0.28 0.27   -0.78   0.13 1.00      779      821
##    z1[6]   -0.51  -0.51 0.28 0.27   -0.96  -0.02 1.00      912      982
##    z1[7]    1.42   1.42 0.33 0.31    0.89   1.95 1.00     1295     1185
##    z1[8]    0.53   0.51 0.46 0.46   -0.21   1.33 1.00      626      809
##    z1[9]    1.33   1.33 0.30 0.29    0.84   1.82 1.00     1358     1248
##    z1[10]  -0.59  -0.59 0.29 0.28   -1.06  -0.11 1.00      982     1190
##    z1[11]  -0.40  -0.41 0.28 0.27   -0.85   0.07 1.00      823      947
##    z1[12]   0.66   0.65 0.24 0.23    0.28   1.06 1.00     1465     1598
##    z1[13]   1.92   1.92 0.49 0.47    1.12   2.74 1.00     1101     1151
##    z1[14]  -0.68  -0.68 0.30 0.30   -1.18  -0.19 1.00     1067     1199
##    z1[15]  -0.16  -0.16 0.29 0.29   -0.62   0.34 1.00      682      747
##    z1[16]   0.42   0.41 0.28 0.28   -0.03   0.89 1.00     1180     1373
##    z1[17]  -0.73  -0.73 0.31 0.30   -1.24  -0.21 1.00     1104     1153
##    z1[18]   0.32   0.31 0.31 0.30   -0.19   0.84 1.00     1092     1485
##    z1[19]   0.29   0.28 0.32 0.31   -0.23   0.81 1.00     1074     1374
##    z1[20]   0.94   0.93 0.23 0.22    0.57   1.32 1.00     1698     1426
##    z1[21]   0.29   0.28 0.31 0.31   -0.22   0.82 1.00     1078     1439
##    z1[22]   0.52   0.51 0.26 0.26    0.10   0.96 1.00     1284     1654
##    z1[23]   1.87   1.86 0.47 0.45    1.11   2.65 1.00     1117     1165
##    z1[24]   0.92   0.92 0.23 0.22    0.55   1.31 1.00     1702     1450
##    z1[25]   1.25   1.25 0.28 0.27    0.79   1.71 1.00     1425     1254
##    z1[26]  -0.75  -0.75 0.31 0.31   -1.27  -0.22 1.00     1117     1111
##    z1[27]  -0.69  -0.69 0.30 0.30   -1.19  -0.19 1.00     1074     1198
##    z1[28]   1.06   1.06 0.24 0.23    0.67   1.47 1.00     1600     1223
##    z1[29]   0.59   0.58 0.25 0.24    0.19   1.01 1.00     1371     1599
##    z1[30]   0.54   0.52 0.26 0.25    0.13   0.97 1.00     1305     1626

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4 5 6 7
##   0 1 2 0 0 0 0 0 0
##   1 1 2 2 0 0 0 0 0
##   2 0 0 3 2 0 0 0 0
##   3 0 0 0 2 0 1 0 0
##   4 0 0 1 1 3 0 0 0
##   5 0 0 0 0 1 4 0 0
##   6 0 0 0 0 1 0 1 0
##   7 0 0 0 0 0 0 0 1
##   8 0 0 0 0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4 5 6 7
##   0 1 0 0 0 0 0 0 0
##   1 1 2 0 0 0 0 0 0
##   2 0 0 2 2 0 0 0 0
##   3 0 0 0 0 0 0 0 0
##   4 0 0 1 0 1 0 0 0
##   5 0 0 0 0 1 0 0 0
##   6 0 0 0 0 0 0 0 0
##   7 0 0 0 0 0 0 0 0
##   8 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4 5 6 7
##   0 0 2 0 0 0 0 0 0
##   1 0 0 2 0 0 0 0 0
##   2 0 0 1 0 0 0 0 0
##   3 0 0 0 2 0 1 0 0
##   4 0 0 0 1 2 0 0 0
##   5 0 0 0 0 0 4 0 0
##   6 0 0 0 0 1 0 1 0
##   7 0 0 0 0 0 0 0 1
##   8 0 0 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 4 5 6 7 8
##   0 1 1 1 0 0 0 0 0 0
##   1 1 2 2 0 0 0 0 0 0
##   2 0 0 1 4 0 0 0 0 0
##   3 0 0 0 2 0 0 1 0 0
##   4 0 0 1 1 1 2 0 0 0
##   5 0 0 0 0 1 2 2 0 0
##   6 0 0 0 0 1 0 1 0 0
##   7 0 0 0 0 0 0 0 0 1
##   8 0 0 0 0 0 0 0 1 0
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 4 5 6 7 8
##   0 1 0 0 0 0 0 0 0 0
##   1 1 2 0 0 0 0 0 0 0
##   2 0 0 1 3 0 0 0 0 0
##   3 0 0 0 0 0 0 0 0 0
##   4 0 0 1 0 0 1 0 0 0
##   5 0 0 0 0 1 0 0 0 0
##   6 0 0 0 0 0 0 0 0 0
##   7 0 0 0 0 0 0 0 0 0
##   8 0 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 4 5 6 7 8
##   0 0 1 1 0 0 0 0 0 0
##   1 0 0 2 0 0 0 0 0 0
##   2 0 0 0 1 0 0 0 0 0
##   3 0 0 0 2 0 0 1 0 0
##   4 0 0 0 1 1 1 0 0 0
##   5 0 0 0 0 0 2 2 0 0
##   6 0 0 0 0 1 0 1 0 0
##   7 0 0 0 0 0 0 0 0 1
##   8 0 0 0 0 0 0 0 1 0
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

gamma regression

objective variable [0,Infinity)

# of samples is N,  
log (a/li)=sum(bj*xji),j=[0,K],i=[1,N]
li=a/exp(sum(bj*xji))
yi~Ga(a,li)  

ex6-4.stan

gamma regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> a;
}
model{
  vector[N] l;
  for(i in 1:N){
    l[i]=a/exp(X[i]*b);
  }
  y~gamma(a,l);
}
n=20
tb=tibble(x1=runif(n,0,2),x2=runif(n,0,2),
          y=rgamma(n,3,3/exp(x1+x2)))

f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-4.stan')

mle=mdl$optimize(data=data)  # it sometimes occur error and stop process
## Initial log joint probability = -839.095 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -48.9895   0.000130897   0.000166455           1           1       20    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -48.99
##      b[1]    -0.09
##      b[2]     1.27
##      b[3]     0.88
##      a        3.59
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -49.92 -49.57 1.59 1.34 -53.02 -48.04 1.00      770     1126
##      b[1]  -0.07  -0.08 0.30 0.30  -0.58   0.41 1.00     1002     1010
##      b[2]   1.28   1.28 0.27 0.25   0.83   1.73 1.00     1082     1105
##      b[3]   0.88   0.87 0.25 0.25   0.49   1.29 1.00     1315     1099
##      a      3.46   3.34 1.12 1.03   1.88   5.50 1.00      932      876

negative binomial regression

The event with probability p do not occur y times till the event occur a times  
(negative binomial distribution has larger variance than poisson distribution)

y~NB(a,p), log(p)=X*b

y~NB(a,l0/(1+l0)) when y~Po(l), l~Ga(a,l0), l0=a/E[l]

ex6-5.stan

negative binomial regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> a;
}
model{
  a~cauchy(0,5);
  y~neg_binomial_2_log(X*b,a);
}
n=20
tb=tibble(x1=runif(n,-1,0),x2=runif(n,-1,0),
          y=rnbinom(n,3,exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -66.829 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -46.4879     0.0203844   0.000489725           1           1       19    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -46.49
##      b[1]    -0.19
##      b[2]    -1.28
##      b[3]    -1.63
##      a        3.20
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -47.32 -46.94 1.52 1.28 -50.26 -45.59 1.01      647      710
##      b[1]  -0.20  -0.16 0.53 0.52  -1.11   0.60 1.01      743      576
##      b[2]  -1.30  -1.29 0.62 0.60  -2.34  -0.32 1.01     1023     1096
##      b[3]  -1.62  -1.60 0.76 0.76  -2.87  -0.44 1.00      910      798
##      a      4.37   3.49 3.40 1.79   1.52   9.38 1.00     1134      826

beta regression

using prior of binomial distribution parameter p[0,1]
y~B(n,p), p~beta(a,b)

m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)

m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*s

ex6-6.stan

beta regression

data{
  int N;
  int K;
  vector[N] p;
  matrix[N,K] X;
}
parameters{
  vector[K] be;
  real<lower=0> s;
}
model{
  vector[N] a;
  vector[N] b;
  for(i in 1:N){
    a[i]=inv_logit(X[i]*be)*s;
    b[i]=(1-inv_logit(X[i]*be))*s;
  }
  p~beta(a,b);
}
n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
          p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$p)
plot(tb$x2,tb$p)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),p=tb$p,X=X)

mdl=cmdstan_model('./ex6-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -23.3907 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22       5.89803   0.000126732    0.00013279           1           1       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__      5.90
##     be[1]    -1.98
##     be[2]     4.66
##     be[3]     1.80
##     s         4.56
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__   5.37   5.73 1.48 1.21  2.38  7.09 1.01      610      867
##     be[1] -2.07  -2.06 0.68 0.68 -3.18 -0.93 1.01      636      804
##     be[2]  4.89   4.86 1.56 1.57  2.34  7.43 1.01      726      631
##     be[3]  1.84   1.87 1.73 1.73 -0.95  4.59 1.00      656      643
##     s      4.64   4.50 1.37 1.30  2.65  7.12 1.00     1644     1434

beta binomial regression

fitting to distribution has larger variance than binomial distribution
y~betaB(n,a,b) when y~B(n,p), p~beta(a,b)

m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)

m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*s

ex6-7.stan

beta binomial regression

data{
  int N;
  int K;
  array[N] int y;
  array[N] int n;
  matrix[N,K] X;
}
parameters{
  vector[K] be;
  real<lower=0> s;
}
model{
  vector[N] a;
  vector[N] b;
  for(i in 1:N){
    a[i]=inv_logit(X[i]*be)*s;
    b[i]=(1-inv_logit(X[i]*be))*s;
  }
  y~beta_binomial(n,a,b);
  s~cauchy(0,5);
}
n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
          p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3),
          n1=floor(runif(n,5,9)),
          y=rbinom(n,n1,p))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),n=tb$n1,y=tb$y,X=X)

mdl=cmdstan_model('./ex6-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -80.2502 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -78.1212   0.000434371    0.00119507      0.5333      0.5333       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__    -78.12
##     be[1]    -1.73
##     be[2]     3.87
##     be[3]     2.41
##     s         3.38
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##     lp__  -78.84 -78.56 1.41 1.31 -81.61 -77.14 1.00      735     1078
##     be[1]  -1.95  -1.92 0.78 0.73  -3.33  -0.75 1.00      919     1085
##     be[2]   4.17   4.20 1.88 1.89   1.05   7.27 1.00     1003      958
##     be[3]   2.86   2.76 2.17 2.13  -0.59   6.64 1.00      853      866
##     s       5.38   4.06 9.63 2.06   1.81  11.55 1.01     1944     1042